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The  objective  of  this  study  is  to  develop  micromechanical  models  for  predicting  the 
stiffness  and  strength  properties  of  textile  composite  materials.  Micromechanical  analysis 
of  textile  composites  is  possible  due  to  the  presence  of  a repeating  unit-cell  or  representative 
volume  element.  The  unit-cell  is  assumed  to  span  the  material  continuously  in  all  three 
dimensions.  On  the  microscale — comparable  to  the  scale  of  the  unit-cell — the  composite  is 
heterogeneous  due  to  the  presence  of  the  reinforcing  yarn  and  the  matrix.  However,  on  the 
macroscale — comparable  to  the  structural  scale — the  composite  is  assumed  to  be 
homogeneous  and  orthotropic.  The  homogeneous  composite  properties  are  then  predicted 
from  the  constituent  material  properties  and  the  yarn  geometry. 

The  highlight  of  this  study  is  a systematic  analysis  of  issues  involved  in  the  finite 
element  based  micromechanics  of  textile  composites.  The  unit-cell  is  discretized  with 
three-dimensional  finite  elements,  and  periodic  boundary  conditions  are  imposed  between 
opposite  end-faces  of  the  unit-cell.  Six  linearly  independent  deformations  are  applied  to  the 
unit-cell.  From  the  forces  acting  on  the  unit-cell  for  each  of  the  six  deformations,  the 


composite  stiffness  matrix  is  obtained.  A similar  procedure  is  followed  to  determine  the 
composite  coefficients  of  thermal  expansion.  The  numerical  procedure  was  tested  by 
applying  it  for  simple  examples,  for  which  the  results  are  known.  The  numerical  results  were 
also  compared  with  existing  models  for  textile  composites.  In  both  cases,  the  results  compare 
very  favorably.  The  finite  element  procedure  is  extended  to  compute  the  thermal  residual 
microstresses  and  to  estimate  the  initial  failure  envelope  for  the  textile  composite. 

An  independent  finite  element  micromechanical  analysis,  analogous  to  the  above,  is 
presented  for  thin  textile  composite  structures  with  few  unit-cells  in  the  thickness  direction. 
In  that  case,  the  composite  is  modeled  as  a homogeneous  plate  to  predict  the  plate  stiffness 
coefficients  and  plate  coefficients  of  thermal  expansion.  It  was  shown  the  plate  properties 
could  not  be  predicted  from  the  corresponding  three-dimensional  properties. 

In  addition,  an  approximate  analytical  procedure  is  presented  to  estimate  the 
composite  thermo-elastic  constants.  The  procedure  called  the  Selective  Averaging  Method 
is  based  on  a judicious  combination  of  stiffness  and  compliance  averaging.  The  method  is 
fast  and  easy  to  implement,  and  suitable  for  parametric  studies. 


VI 


CHAPTER  1 
INTRODUCTION 

l.I  Background 

The  increasing  demand  for  lightweight  yet  strong  and  stiff  structures  has  lead  to  the 
development  of  advanced  fiber  reinforced  composites.  These  materials  are  not  only  used  in 
the  aerospace  industry  but  also  in  a variety  of  commercial  applications  like  automobile, 
marine  and  biomedical  applications.  Traditionally  fibrous  composites  are  manufactured  by 
laminating  several  layers  of  unidirectional  fiber  tapes  (Fig.  1)  pre-impregnated  with  matrix 
material.  The  effective  properties  of  the  composite  can  be  controlled  by  changing  several 
parameters  like  the  fiber  orientation  in  a layer,  stacking  sequence,  fiber  and  matrix  material 
properties  and  fiber  volume  fraction.  However,  the  manufacture  of  fibrous  laminated 
composites  from  prepregs  is  labor  intensive.  Besides,  laminated  composites  lack 
through-the-thickness  reinforcement,  and  hence  have  poor  interlaminar  strength  and 
fracture  toughness. 


Figure  1.  Stacking  of  layers  (plies)  for  fibrous  laminated  composite. 
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Figure  2.  Examples  of  woven,  braided  and  knitted  textile  preforms.  (Source: 
Chou  and  Ko,  1989) 


Recent  developments  in  textile  manufacturing  processes  show  some  promise  in 
overcoming  the  above  limitations.  Textile  processes  such  as  weaving,  braiding  and  knitting 
can  turn  large  volumes  of  yam  into  dry  preforms  at  a faster  rate,  thus  reducing  costs  and  cycle 
time.  The  dry  preforms  (Fig.  2)  are  impregnated  with  appropriate  matrix  material  and  cured 
in  a mold  using  processes  such  as  Resin  Transfer  Molding  (RTM).  Two-dimensional  woven 
and  braided  mats  offer  increased  through  the  thickness  properties  due  to  yarn  interlacing. 
The  mats  may  be  stitched  using  Kevlar  or  glass  threads  to  provide  additional  reinforcement 
along  the  thickness  direction  (Sharma  and  Sankar,  1995).  Three-dimensional  woven  and 
braided  composites  provide  multidirectional  reinforcement,  thus  directly  enhancing  the 
strength  and  stiffness  in  the  thickness  direction.  Unlike  laminated  structures 
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three-dimensional  composites  do  not  possess  the  weak  plane  of  delamination,  thus  giving 
increased  impact  resistance  and  fracture  toughness.  Textile  manufacturing  processes  in 
conjunction  with  resin  transfer  molding  are  also  suitable  for  the  production  of  intricate 
structural  forms  at  a reduced  cycle  time.  This  allows  complex  shaped  structures  to  be 
fabricated  as  an  integral  unit,  thus  eliminating  the  use  of  Joints  and  fasteners. 

With  the  advancements  in  aforementioned  technologies  there  is  a need  to  develop 
scientific  methods  of  predicting  the  performance  of  the  composites  made  using  the  above 
processes.  There  are  numerous  variables  involved  in  textile  processes  besides  the  choice  of 
the  fiber  and  matrix  materials.  This,  for  example,  includes  (a)  the  number  of  filaments  in  the 
yarn  specified  by  the  yarn  linear  density  and  (b)  the  yarn  architecture  (description  of  the  yarn 
geometry)  determined  by  the  type  of  weaving  or  braiding  processes.  Thus,  there  is  a need 
for  analytical/numerical  models  to  study  the  effect  of  these  variables  on  the  textile  composite 
behavior. 

Ideally  a structural  engineer  would  like  to  model  textile  composites  as  a 
homogeneous  anisotropic  material — preferably  orthotropic — so  that  the  structural 
computations  can  be  simplified,  and  also  the  existing  computer  codes  can  be  used  in  the 
design.  This  would  require  the  prediction  of  the  effective  (macroscopic)  properties  of  the 
composites  from  the  constituent  material  (microscopic)  characteristics  such  as  yam  and 
matrix  properties,  yam-matrix  interface  characteristics  and  the  yarn  architecture  (Fig.  3). 
This  is  possible  if  we  assume  that  there  is  a representative  volume  element  (RVE)  or  an 
unit-cell  that  repeats  its  self  throughout  the  volume  of  the  composite,  which  is  tme  in  the  case 
of  textile  composites.  The  unit-cell  can  be  considered  as  the  smallest  possible  building  block 
for  the  textile  composite,  such  that  the  composite  can  be  created  by  assembling  the  unit-cell 
in  all  three  dimensions  (Fig.  4).  The  prediction  of  the  effective  macroscopic  properties  from 
the  constituent  material  characteristics  is  one  of  the  aspects  of  the  science  known  as 
micromechanics.  The  effective  properties  include  thermo-mechanical  properties  like 
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stiffness,  strength  and  coefficients  of  thermal  expansion  as  well  as  thermal  conductivities, 
electromagnetic  and  other  transport  properties. 


Figure  3.  Essence  of  micromechanics. 

(a)  composite  (microscale)  showing  plain-weave  pattern 

(b)  homogeneous  and  orthotropic  composite  (macroscale). 
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Most  of  the  early  work  to  determine  the  properties  of  heterogeneous  materials  was 
restricted  to  particulate  (spherical)  inclusions  which  were  assumed  to  be  isotropic.  For 
example,  Hashin  (1962)  derived  expressions  for  the  bounds  for  the  elastic  moduli  of 
heterogeneous  materials  using  variational  theory.  A review  of  analytical  methods  for 
predicting  the  effective  properties  of  particulate  composites  was  presented  by  Christensen 
(1990).  However,  the  complex  geometry  of  the  textile  preforms  makes  such  precise 
analytical  modeling  difficult. 

The  current  literature  dealing  with  micromechanical  analyses  for  textile  composites 
can  be  broadly  classified  into  three  categories:  mechanics  of  materials  type  models,  energy 
based  approach,  and  finite  element  analysis  of  the  unit-cell.  All  of  the  above  models 
recognize  that  there  is  a unit-cell  in  the  composite  material,  and  they  attempt  to  model  the 
material  as  a homogeneous  orthotropic  material.  In  the  mechanics  of  materials  type  models, 
the  yarns  are  approximated  as  simple  structural  elements,  e.g. , beams,  plates,  laminates  etc., 
and  their  deformation  behavior  is  assumed  to  be  governed  by  the  corresponding  structural 
constitutive  relations.  The  kinematics  is  also  simplified  to  a great  extent,  and  a relation 
between  the  overall  deformation  of  the  unit-cell  and  the  average  forces  are  derived.  The 
energy  approach  is  similar  to  the  previous  one,  except  that  the  strain  energy  in  the  unit-cell 
is  evaluated  based  on  some  assumed  displacement  field,  which  is  usually  an 
oversimplification  of  exact  displacements.  The  elastic  constants  are  derived  by  equating  the 
strain  energy  in  the  approximate  model  to  that  in  an  idealized  homogeneous  composite. 
Most  energy  based  approaches  provide  bounds  for  the  homogeneous  properties,  and  can  be 
used  as  a check  for  experimental  observations  or  other  theoretical  models.  The  third  method 
is  the  rigorous  micromechanical  analysis  of  the  unit-cell,  which  often  requires  the  use  of 
numerical  methods  such  as  the  finite  element  method,  and  also  uses  the  exact 
three-dimensional  constitutive  relations  for  the  yarn  and  the  matrix  material. 

Ishikawa  and  Chou  (1982a,  1982b,  1983a,  1983b,  1983c)  proposed  three  analytical 
models  for  thermo-elastic  properties  of  woven  fabric  composites — the  mosaic  model,  the 
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fiber  undulation  model  and  bridging  model.  The  mosaic  model  treats  the  composite  as  an 
assembly  of  crossply  laminates,  and  then  uses  lamination  theory  to  predict  the  composite 
properties.  The  fiber  undulation  model  in  addition  takes  into  account  the  effect  of  yam 
continuity  in  the  loading  direction.  The  bridging  model  was  specifically  developed  to 
estimate  the  elastic  behavior  of  satin  weaves.  This  model  simulates  the  effect  of  load 
distribution  and  load  transfer  between  the  yarns.  The  above  three  models  were  developed 
for  two-dimensional  fabric  composites.  These  models  were  then  extended  by  Yang  et  al. 
(1986)  in  the  fiber  inclination  model  to  predict  the  elastic  properties  of  three-dimensional 
textile  composites.  The  fiber  inclination  model  was  used  to  determine  the  elastic  properties 
of  three-dimensional  angle-interlock  composites  and  braided  composites  by  Whitney  and 
Chou  (1989)  and  Ma  et  al.  (1986)  respectively.  The  results  were  compared  with  test  data, 
and  parametric  analyses  were  performed  to  study  the  effect  of  changing  geometric 
parameters.  Crane  and  Camponeschi  ( 1 986)  presented  an  analytical  model  based  on  classical 
lamination  theory  to  predict  the  extensional  stiffnesses  for  multi-directional  braided 
composites.  Naik  (1994)  used  a stiffness  averaging  technique  to  predict  woven  and  braided 
composite  properties  (explained  in  Chapter  3). 

Numerical  modeling  of  the  unit-cell  is  popular  due  to  its  ability  to  capture  the  effects 
of  complicated  yarn  architectures.  For  instance,  Yoshino  and  Ohtsuka  (1982)  performed  a 
two-dimensional  finite  element  analysis  using  plane  strain  elements  to  predict  the  stress 
distribution  within  a plain-weave  fabric.  Dasgupta  et  al.  (1990)  and  Whitcomb  (1991) 
analyzed  the  unit-cell  of  a plain-weave  composite  using  three-dimensional  finite  elements 
to  determine  the  effect  of  the  yarn  geometry  and  yarn  volume  fraction  on  the  composite 
thermo-elastic  constants.  Foye  (1993)  used  inhomogeneous  elements  called  replacement 
elements  to  model  the  unit-cell.  His  model  can  be  used  to  predict  both  composite  stiffness 
and  strength  properties.  Cox  et  al.  (1994)  presented  a three-dimensional  finite  element  model 
using  two  types  of  elements.  The  yams  are  modeled  as  two-node  line  elements  and  the  rest 
of  the  medium  as  eight-node  solid  elements.  The  model  was  used  to  predict  failure 
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mechanisms  in  angle  and  orthogonal  interlock  woven  composites.  Raju  et  al.  (1990) 
compiled  a review  of  available  analytical  and  numerical  methods  for  modeling  textile 
composites. 


1 .2  Scope  of  Study 

The  analytical  models  explained  above  are  approximate — they  make  use  of 
assumptions  similar  to  that  in  classical  plate  theory,  and  the  yarn  geometry  is  greatly 
simplified.  Further,  these  models  do  not  consider  periodic  boundary  conditions  (see  Section 
2.1.1)  for  continuity  of  displacements  and  tractions  on  opposite  faces  of  the  unit-cell.  Most 
of  the  existing  numerical  models  are  for  specific  yarn  geometries  such  as  simple  weave 
architectures.  Therefore  there  is  a need  for  more  general  models  to  predict  the  mechanical 
properties  for  a textile  composite,  particularly  for  composites  with  complex  woven  and 
braided  geometries.  The  present  study  illustrates  two  such  methods.  The  first  method 
involves  rigorous  finite  element  analysis  of  the  unit-cell,  which  imposes  exact  displacement 
and  traction  boundary  conditions  on  the  end-faces  of  the  unit-cell.  The  second  is  an 
approximate  analysis  in  which  the  unit-cell  is  discretized  into  a number  of  elements.  The 
effective  (macroscale)  stiffness  for  the  composite  is  computed  by  a combination  of  stiffness 
and  compliance  averaging  of  the  element  stiffness  coefficients.  The  advantage  of  the 
proposed  methods  is  their  generality — they  are  applicable  to  any  textile  geometry  provided 
that  there  exists  a repeatable  pattern  (which  defines  the  unit-cell).  Both  the  methods  can 
compute  the  effective  elastic  constants  and  effective  coefficients  of  thermal  expansion  for 
the  homogenous  composite. 

The  analysis  for  the  elastic  constants  of  the  textile  composite  continuum  assumes  that 
the  unit-cell  is  repeatable  in  all  three  dimensions.  This  assumption  will  not  be  valid  for  thin 
textile  composite  structures,  where  there  will  be  a only  a few  unit-cells  in  the  thickness 
direction.  For  such  applications,  an  alternate  micromechanical  analysis  is  proposed.  The 
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unit-cell  is  modeled  as  a plate  (Fig.  5)  and  the  corresponding  plate  stiffness  coefficients  and 
plate  coefficients  of  thermal  expansion  are  predicted. 


Chapter  2 explains  the  finite  element  procedure  to  compute  the  macroscale 
thermo-elastic  constants.  The  procedure  is  detailed  for  a thick  textile  composite  in  which 
there  are  sufficient  number  of  unit-cells  in  the  thickness  direction.  The  procedure  is  then 
extended  to  a thin  composite.  The  thin  composite  is  first  modeled  as  a beam  and  a 
two-dimensional  analysis  is  presented  to  compute  the  beam  thermo-elastic  coefficients.  A 
similar  three-dimensional  analysis  is  also  illustrated  in  which  the  thin  composite  is  modeled 
as  a plate.  An  analytical  procedure  called  the  Selective  Averaging  Method  (SAM)  is 
presented  in  Chapter  3 to  estimate  the  macroscale  thermo-elastic  constants.  This  procedure 
utilizes  a judicious  combination  of  stiffness  and  compliance  averaging  to  compute  the 
composite  properties.  Again  for  a thin  textile  composite  the  plate  stiffness  coefficients  and 
plate  coefficients  of  thermal  expansion  may  be  computed. 
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In  Chapter  4 the  microstresses  computed  using  the  finite  element  analysis,  are  used 
to  construct  an  initial  failure  envelope  for  a textile  composite  modeled  as  homogeneous 
continuum.  An  independent  procedure  using  plate  theory  is  presented  to  predict  the  failure 
envelope  for  a thin  textile  composite.  Since  the  composite  is  generally  fabricated  at  a 
temperature  greater  than  room  temperature,  thermal  residual  microstresses  are  developed 
due  to  the  mismatch  in  the  coefficients  of  thermal  expansion  for  the  constituent  materials. 
A method  to  predict  these  thermal  microstresses  is  illustrated.  Chapter  5 deals  with  the  issues 
pertaining  to  finite  element  modeling  of  the  unit-cell.  Due  to  the  complex  yam  architectures, 
traditional  finite  element  modeling  becomes  exceedingly  difficult.  Various  alternative 
means  of  finite  element  modeling  the  unit-cell  are  discussed  with  suitable  two-dimensional 
examples.  Finally  Chapter  6 summarizes  the  work  done  in  this  study.  The  limitations  of  the 
work  are  addressed  and  suggestions  are  given  for  future  work  in  this  area. 


CHAPTER  2 

FINITE  ELEMENT  MODELS  FOR  THERMO-ELASTIC  CONSTANTS 


In  this  chapter,  we  demonstrate  finite  element  based  micromechanical  models  to 
predict  the  effective  stiffness  properties  and  effective  coefficients  of  thermal  expansion 
(CTE’s)  for  a textile  composite.  The  macroscale  properties  of  the  composite  are  determined 
at  a scale  much  larger  than  the  dimensions  of  the  unit-cell,  but  comparable  to  the  dimensions 
of  the  structural  component.  The  average  stresses  at  a point  at  the  structural  scale  will  be 
called  the  macroscale  stresses  or  macrostresses.  The  actual  stresses  at  a point  at  the 
continuum  level  will  be  called  the  microscale  stresses  or  microstresses.  To  distinguish  the 
macroscale  deformations  and  stresses  from  their  microscale  counterparts — a superscript 
“M”  will  be  used  to  denote  the  macroscale  deformations  and  stresses. 

2. 1 Unit-Cell  Analysis  for  Three-Dimensional  Elastic  Constants 

In  this  section  a procedure  to  determine  the  three-dimensional  elastic  constants  for 
a textile  composite  material  is  described.  Consider  a rectangular  hexahedron  as  the  unit-cell 
of  the  three-dimensional  textile  composite.  The  edges  of  the  unit-cell  are  assumed  to  be 
parallel  to  the  coordinate  axes  xj,  X2  andxj,  with  unit-cells  repeating  in  all  three  directions. 
The  length  of  the  unit-cell  in  the  X[  direction  is  defined  as  L, . On  the  macroscale  the  composite 
is  assumed  to  be  homogeneous  and  orthotropic  and  the  composite  behavior  is  characterized 
by  the  following  constitutive  relation  : 
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"^23 
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*'11 
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^31 

yfi 

k.  .. 
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where  {o^]  and  {e^}  are  the  macrostresses  and  macrostrains,  respectively;  {a^}  and  [C] 
are  the  macroscale  CTE’s  and  orthotropic  elasticity  matrix  to  be  determined;  AT^  is  a 
uniform  temperature  difference  throughout  the  unit-cell.  The  macroscale  elasticity  matrix 
is  of  the  form: 


Cii  C^2 

^13 

Ci4 

^15 

C,6 

^22 

^23 

C24 

^25 

^26 

^33 

^4 

^35 

^36 

symm. 

C44 

Qs 

^55 

^56 

^66 


(2) 


The  unit-cell  analysis  assumes  that  the  material  is  homogeneous  and  orthotropic  on  the 
macroscale.  The  material,  therefore,  is  subject  to  a uniform  state  of  strain  in  the  macroscopic 
sense.  The  macrostresses  required  to  create  such  a state  of  strain  are  computed  from  the  finite 
element  model  of  the  unit-cell.  In  the  microscale,  all  unit-cells  have  identical  stress  and  strain 
fields.  Continuity  of  microstresses  across  the  unit-cell  then  requires  that  tractions  be  equal 
and  opposite  at  corresponding  points  on  opposite  faces  of  the  unit-cell.  Since  the 
displacement  gradients  are  constant  for  a homogeneous  deformation,  the  displacements  at 
corresponding  points  on  opposite  faces  of  the  unit-cell  differ  only  by  a constant. 
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2.1.1  Periodic  Boundary  Conditions 

The  periodic  boundary  conditions  (EC’s)  consist  of  the  periodic  displacement 
boundary  conditions  which  ensure  the  compatibility  of  displacements  on  opposite  faces  of 
the  unit-cell,  and  the  periodic  traction  boundary  conditions  to  enforce  the  continuity  of 
stresses.  A macroscopically  homogeneous  deformation  can  be  represented  as 

uf  = H^jXj  ij  =1,2,3  (3) 

where  Hy  are  the  displacement  gradients.  Then  the  periodic  displacement  boundary 
conditions  to  be  imposed  on  the  faces  Xi=0  and  Xi=L{  are 

m-(LjX2_a:3)  — m,(0,j:2,^3)  = 

Uiix^  L2  X2)  - M,(xi  0,JC3)  = //,-2^2  (4) 

-^2^3)  ~ u■(x^X20)  = H-^L^ 

The  traction  boundary  conditions  to  be  imposed  on  the  faces  jc,=0  and  Xi=L\  are 

FIL^  X2,X2)  = - F,(0,X2,JC3) 

F,(xi  L2,jc3)  = - F,(.x,  0,^3)  (5) 

F,(ac,  X2L3)  = - F.(Xj  X2,0) 

The  periodic  displacement  and  traction  EC’s  (for  a two-dimensional  problem)  to  be  imposed 
for  the  deformation  given  by  £77  =i/L;  and  £22=712=0  are  shown  in  Fig.  6. 

The  above  periodic  EC’s  are  imposed  in  the  finite  element  model  either  by  using 
multi-point  constraint  elements  or  by  incorporating  transformation  equations  to  eliminate 
the  constrained  displacements  (Cook  et  al.,  1989).  These  two  methods  to  impose  the  periodic 
EC’s  are  discussed  in  Section  5.1.3.  Eoth  the  methods  require  a finite  element  model  with 
corresponding  nodes  on  opposite  faces  of  the  unit-cell. 
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2.1.2  Determination  of  Three-Dimensional  Elastic  Constants  and  CTE’s 

The  unit-cell  is  discretized  with  three-dimensional  finite  elements  such  that  opposite 
faces  of  the  unit-cell  have  identical  nodes.  Periodic  displacement  and  traction  boundary 
conditions  are  enforced  between  the  corresponding  nodes.  The  periodic  displacement  EC’s 
are  imposed  such  that  only  one  of  the  macrostrain  components  is  non-zero;  and  the  uniform 
temperature  difference  A is  set  to  zero.  Then,  the  difference  in  displacements  between 
corresponding  points  on  opposite  faces  of  the  unit-cell  will  be  equal  to  that  in  a homogenous 
continuum  subject  to  the  same  deformation.  The  average  stresses  (macrostresses)  required 
to  create  such  a deformation  are  obtained  from  the  finite  element  results.  Substituting  the 
macrostresses  and  macrostrains  in  the  composite  constitutive  relation  (Eqn.  1)  the  stiffness 
coefficients  in  the  column  corresponding  to  the  non-zero  strain  can  be  evaluated.  This 
procedure  is  repeated  for  other  macrostrain  components  (keeping  the  temperature  difference 
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zero)  to  obtain  the  entire  stiffness  matrix  [C].  The  orthotropic  elastic  constants  of  the 
composite  material  can  be  easily  determined  by  inverting  the  stiffness  matrix,  and 
comparing  the  compliance  coefficients  with  that  of  an  orthotropic  material. 

To  compute  the  six  CTE’s,  a finite  temperature  change  To  is  applied  to  all  the 
elements  in  the  unit-cell;  and  periodic  displacement  EC’s  are  imposed  such  that  all  the 
macrostrain  components  are  zero.  Then  the  composite  constitutive  relation  Eqn.  (1)  will 
reduce  to 

[a^]  = - To  (6) 

The  macrostresses  for  such  a deformation  are  computed  as  described  below.  Knowing  the 
stiffness  matrix  [C],  the  composite  CTE’s  are  found  as 

|a")  = -:p[C]-’(a^)  (7) 

^ O 


Table  1 presents  the  non-zero  displacement  EC’s  imposed  on  the  unit-cell  to  obtain  [C]  and 
the  CTE’s 

Table  1.  Non-zero  displacement  EC’s  to  obtain  3-D  elastic  constants  and  CTE’s. 


stiffness  coefficients  to  be  obtained 

non-zero  displacement  EC’s  | 

first  column  of  [C]  (£77^  = 1) 

U](Lj,  X2,  X3)  -ui(0,  X2,  xj)  = L] 

second  column  of  [C]  (£22^  = 1 ) 

U2(X],  L2,  X3)  - U2(xj,  0,  X3)  = L2 

third  column  of  [C]  (£jj^  = 1 ) 

U3(x],  X2,  L3)  - U3(x],  X2,  0)  = L3 

fourth  column  of  [C]  (y2S^  = 

U2(X],  X2,  L3)  - U2(X],  X2,  0)  = L3/2 
U3(X],  L2,  X3)  - U3(X],  0,  X3)  = L2/2 

fifth  column  of  [C]  (y/s^  = 1) 

uj(x],  X2,  L3)  - ui(x],  X2,  0)  = L3/2 
U3(L],  X2,  X3)  - U3(0,  X2,  X3)  = Lj/2 

sixth  column  of  [C]  (772^  = 1) 

U](xj,  L2,  X3)  - U](xi,  0,  X3)  = L2/2 
U2(L],  X2,  X3)  - U2(0,  X2,  X3)  = Lj/2 

CTE’s  (zir^=  1) 

AT=  1 
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The  macrostresses  for  a given  deformation  state  can  be  found  by  the  following  two 
methods.  In  the  first  method,  the  macrostresses  are  obtained  by  averaging  the  nodal  forces 
on  each  face  of  the  unit  cell.  For  example,  the  macroscale  can  be  obtained  as 


= 

^11 


L2T3 


(8) 


where  is  the  nodal  force  in  the  xj  direction  at  the  nth  node,  and  ^ denotes  summation 

n 

over  all  nodes  on  the  face  xj  =Lj . Alternatively,  the  macrostresses  can  be  computed  by 
volume-averaging  the  corresponding  microstress  component  in  the  unit-cell.  Then  the 
macrostress  component  cr^j  is  obtained  as 


a 


M 

11 


1 


Oi]{x,y,z)dV 


(9) 


V 

where  V is  the  volume  of  the  unit-cell.  The  microstresses  are  computed  at  the  quadrature 
points,  and  numerically  integrated  over  the  volume  in  each  element  of  the  unit-cell. 


2.1.3  Results  and  Discussion  for  3-D  Elastic  Constants  and  CTE’s 


The  above  procedure  was  demonstrated  for  the  following  material  systems: 
Example  1 . isotropic  material 

Example  2.  bimaterial  medium — both  materials  are  assumed  isotropic  (Fig.  7) 

Example  3 . unidirectional  composite  with  identical  Poisson  ratios  for  fiber  and  matrix — 

fiber  and  matrix  materials  are  isotropic  (Fig.  8) 

Example  4.  unidirectional  composite  with  different  Poisson  ratios  for  fiber  and  matrix — 
fiber  and  matrix  materials  are  isotropic 

Example  5.  plain-weave  textile  composite  (Fig.  9) — yarn  geometry  and  properties 
obtained  from  Dasgupta  et  al.  (1990) 

Example  6.  plain-weave  textile  composite — yam  geometry  and  properties  obtained  from 

Naik(1994) 
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Example  7.  5-hamess  satin  weave  (Fig.  10) — yam  geometry  and  properties  obtained 
from  Naik  (1994) 

For  the  textile  composite  examples  (examples  5-7),  the  yarn  is  assumed  to  be  transversely 
isotropic  and  the  matrix  material  is  assumed  isotropic.  The  constituent  material  properties 
for  the  examples  are  listed  in  Table  2. 


Table  2.  Properties  of  constituent  materials  for  examples  1-7. 


Example  1 

unit-cell  size: 

£ = 70  GPa,  v=0.3,  a=10xl0-^/<>C 
0.500x0.500x0.256  mm 

Example  2 

layer  1 (E-glass): 
layer  2 (epoxy): 
unit-cell  size: 

Ej=70  GPa,  VI  =0.200,  aj=  5x10-^ /<>C,  Vi=0.5 
E2=3.50  GPa,  V2  =0.350,  02  = 60x10-^ /<>€,  V2=0.5 
0.500x0.500x0.256  mm 

Example  3 

fiber: 
matrix: 
unit-cell  size: 

Ef=100  GPa,  Vf=0.300,  a/=  10x10-^ /<>C,  Vf=0.6 
Em  = 10  GPa,  Vm=0.300,  am=  100x10-^ /°C 
10x10x10  fim 

Example  4 

fiber  (E-glass): 
matrix  (epoxy): 
unit-cell  size: 

Ef=70  GPa,  Vf=0.200,  af=  5xlO-^/‘>C,  Vf=0.6 
Em=3.50  GPa,  v„=0.350,  0^=  60x10-^ /^C 
10x10x10  fim 

Example  5 

yarn  (glass-epoxy): 

Ei=58.6l  GPa,  Et=14.49  GPa,  Glt=5.38  GPa,  vlt=0.250 
vrr=0.247,  01=6.1 5x10-^ /‘>C,  ot=22.64x10-^/<’C,  Vf=0.26 

matrix  (epoxy): 

E=3.45GPa,v=0.37,  a=69xI(P^/^C 
unit-cell  size:  1.680x1.680x0.228  mm 

Examples  6,  7 

yarn  (graphite-epoxy): 

Ei= 144.80  GPa,  Et=11.73  GPa,  Glt=5.52  GPa,  Vlt=0.230 
vrr=0.300,  oi=-0.324xl0-^/‘>C,  ot=1  4.00x1 0-^/<>C,  Vf=0.64 

matrix  (epoxy): 

E=3.45  GPa,  v=0.35,  a=40xl0-^/"C 

unit-cell  size:  2.822x2.822x0.2557  mm  (Example  6) 

7.055x7.055x0.2557  mm  (Example  7) 

Vi  stands  for  the  volume  fraction  of  the  constituent  material 
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Figure  7.  Bimaterial  example. 

(a)  bimaterial  medium;  (b)  unit-cell. 


(a)  (b) 


Figure  8.  Unidirectional  composite  example, 
(a)  composite;  (b)  unit-cell. 
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Figure  10.  Yarn  pattern  in  a 5-harness  satin  weave  preform,  (unit-cell  boundary 
in  dotted  lines) 
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A 3-D  finite  element  code  called  fiTEx-XQ  (pronounced  as  microtech)  was  written 
and  implemented  for  the  seven  examples  to  compute  the  homogeneous  elastic  constants  and 
CTE’s.  The  unit-cell  was  assumed  to  be  a rectangular  hexahedron  with  edges  along  the  x-, 
y-  and  z-  axes.  The  unit-cell  was  divided  into  uniform-sized  eight-node  brick  elements  as 
shown  in  Figure  11.  The  number  of  divisions  along  the  x-,  y-  and  z-  axes  are  specified  by 
the  user.  The  elements  were  in  general,  inhomogeneous  elements,  i.e. , consisted  of  more  than 
one  constituent  material.  While  computing  the  element  stiffness  matrix,  the  material 
property  at  the  Gauss  integration  points  was  determined,  and  the  corresponding  elasticity 
matrix  was  used  to  perform  the  integration  over  the  element  volume  (see  Section  5.1 .2).  Thus 
the  element  stiffness  matrix  was  obtained  as  the  averaged  stiffness  of  the  different  materials 
in  the  element.  Periodic  displacement  and  traction  boundary  conditions  were  imposed 
between  corresponding  nodes  on  opposite  faces  of  the  unit-cell.  The  numerical  procedure 
to  compute  the  element  stiffness  matrix  and  to  impose  periodic  EC’s  are  discussed  in  Chapter 
5.  The  skyline  solver  (Bathe,  1982)  was  used  to  solve  for  the  nodal  displacements.  The 
computed  elastic  constants  for  the  seven  examples  are  listed  in  Tables  3-5. 


(b) 


LJy 


Figure  11.  3-D  finite  element  mesh  to  compute  elastic  constants, 
(a)  unit-cell;  (b)  eight-node  brick  element. 
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The  ^TEx-\Q  code  was  first  checked  by  computing  the  elastic  constants  for  an 
isotropic  medium  (Table  3).  Then  the  code  was  implemented  for  a bimaterial  medium.  The 
bimaterial  medium  consisted  of  two  different  layers  of  equal  thickness  in  the  xy-plane 
alternating  in  the  z-direction  (Fig.  7).  The  effective  Young’s  moduli,  Poisson  ratios  and 
CTE’s  of  the  bimaterial  medium  were  derived  exactly,  as  described  below.  The  constitutive 
relation  (considering  only  the  normal  stresses)  for  each  layer  is  defined  as: 


0‘xx 

^11 

C' 

''"12 

C' 

'-13 

C 'S 

p} 

^XX 

Oyy 

> — 

r' 

*^12 
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'-'22 
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'-23 

< 

tyy 

i = 1,2 

o‘zz 

J 

0 

'-'13 

C' 

''"23 

*-33 

^ZZ 

(10) 


where  the  superscript  refers  to  the  layer  number.  To  derive  Cj],Cj2  and  C75  for  the  bimaterial 
medium,  was  assumed  as  one;  and  Eyy^  and  were  assumed  to  be  zero.  The 
assumption  of  Szz^=0  and  the  fact  that  the  layers  are  of  equal  thickness  imply  that  = 
-Ezz^.  The  following  constraints  were,  in  addition,  imposed  across  the  bimaterial  interface 


II 

_ pM 

^xx 

tyy 

II 

_ pM 

tyy 

(11) 

= olz 

= 

^zz 

From  the  above  interfacial  constraints  and  Eqn.  (10),  the  stresses  in  each  layer  were 
computed.  The  stresses  in  the  layers  were  volume-averaged  to  yield  the  corresponding 
macrostresses,  i.e.,  Oxx^,  Oyy^  and  Since  was  equal  to  unity,  the  computed 
macrostresses  were  identical  to  the  stiffness  coefficients  in  the  first  eolumn,  namely,  C/y,  C72 
and  Cyj.  A similar  procedure  was  followed  to  find  the  remaining  stiffness  coefficients  and 
CTE’s  for  the  bimaterial  medium.  The  inplane  shear  modulus  of  the  bimaterial  medium  was 
computed  as  Gxy  = (Gy  + G2 )/2  knowing  that  the  shear  strain  was  uniform  in  both  layers.  The 

2G  jG2 

isostress  assumption  was  used  to  derive  the  transverse  shear  modulus  as  Gxz  = -q — q • 
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It  was  found  that  nTEx~\Q  results  were  identical  to  the  elasticity  results  for  the  bimaterial 
medium  (Table  3). 


Table  3.  Continuum  properties  for  examples  1 and  2 using  finite  elements. 


Ex,  Ey 

(GPa) 

Ez 

(GPa) 

Gxz, 

^yz 

(GPa) 

Gjcy 

(GPa) 

^yz 

X 

10-6/oC 

X 

io-6/°c 

Example  1 
(isotropic 

liTEx-\Q 

(FEA) 

10 

10 

3.85 

3.85 

0.300 

0.300 

10 

10 

medium) 

exact 

solution 

10 

10 

3.85 

3.85 

0.300 

0.300 

10 

10 

Example  2 
(bimaterial 

fiTEx-\Q 

(FEA) 

36.79 

9.79 

2.48 

15.23 

0.312 

0.208 

8.19 

59.60 

medium) 

exact 

solution 

36.79 

9.79 

2.48 

15.23 

0.312 

0.208 

8.19 

59.60 

Table  4.  Continuum  properties  for  examples  3 and  4 using  finite  elements. 


El 

(GPa) 

Ej 

(GPa) 

Glt 

(GPa) 

Gjj 

(GPa) 

vlt 

Vtt 

a-L 

xlO-^ 

l°C 

df 

xlO-6/°C 

Example  3 

fiTEx-lO 

(FEA) 

63.55 

36.48 

12.93 

9.94 

0.300 

0.232 

15.74 

40.79 

(uni 

directional 

composite) 

rule  of 
mixt./ 
Halpin- 
Tsai  eqns. 

64 

34.55 

11.26 

13.29 

0.300 

0.300 

15.63 

55.11 

Example  4 

fzTEx-lO 

(FEA) 

43.12 

18.15 

5.59 

3.92 

0.242 

0.222 

7.40 

25.44 

(uni 

directional 

composite) 

rule  of 
mixt./ 
Halpin- 
Tsai  eqns. 

43.40 

14.79 

4.45 

5.91 

0.260 

0.252 

6.77 

34.24 

Table  4 presents  the  elastic  constants  and  CTE’s  for  the  two  unidirectional  composite 
examples.  The  unidirectional  composite  unit-cell  is  shown  in  Fig.  8.  The  unidirectional 
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composite  properties  were  compared  with  available  analytical  solutions.  The  rule  of 
mixtures  formulae  were  used  to  predict  Ei  and  Vu  ; the  Halpin-Tsai  equations  (Halpin  and 
Tsai,  1969)  for  Ej,  Gu  and  vjt  and  Schapery’s  expressions  (Agarwal  and  Broutman,  1990) 
for  the  thermal  coefficients  and  aj.  To  compare  Gjt,  transverse  isotropy  was  assumed  in 


the  T-T  plane  and  computed  as 


2(1  + ^7t) 


However,  the  finite  element  results  show  that 


the  assumption  of  isotropy  in  the  transverse  plane  is  not  valid.  The  assumption  would  give 
Gjt  as  14.81  GPa  and  7.43  GPa  for  examples  3 and  4 respectively,  thus  grossly 
over-estimating  the  finite  element  results.  The  relations  for  Ej,  Vu  and  aj  are  exact  when 
the  poisson  ratios  are  identical  for  the  fiber  and  the  matrix. 


Table  5.  Continuum  properties  for  examples  5,  6 and  7 using  finite  elements. 


Ex,  Ey 
(GPa) 

Ez 

(GPa) 

Gxz, 

Gyz 

(GPa) 

Gxy 

(GPa) 

Vxz^ 

'^yz 

Vxy 

o-x'^^o-y^ 

xl0-6/° 

C 

10-6/oc 

Example  5 

(plain- 

weave) 

^lTEx-\0 

(FEA) 

11.81 

6.14 

1.84 

2.15 

0.408 

0.181 

28.36 

79.57 

Dasgupta 

results 

14.38 

6.25 

1.94 

3.94 

0.463 

0.167 

22.50 

86.00 

Example  6 

(plain- 

weave) 

nTEx-\Q 

(FEA) 

53.61 

10.88 

4.41 

4.72 

0.365 

0.128 

1.56 

22.1  \ 

TEXCAD 

64.38 

11.49 

5.64 

4.87 

0.396 

0.027 

1.33 

20.71 

Example  7 

(5-harness 

weave) 

fiTEx-\0 

(FEA) 

64.51 

11.33 

4.45 

4.85 

0.329 

0.047 

1.55 

22.03 

TEXCAD 

66.33 

11.51 

4.93 

4.89 

0.342 

0.034 

1.46 

21.24 

Figures  9 and  10  illustrate  the  weave  patterns  for  the  plain-weave  (examples  5 and 
6)  and  the  5-harness  satin  weave  (example  7)  respectively.  The  properties  for  example  5 were 
compared  with  Dasgupta’s  results  for  an  overall  fiber  volume  fraction  {Vf)  of  0.26.  The  yam 
properties  were  not  specified  in  Dasgupta  et  al.  (1990).  So  the  rule  of  mixtures  and 
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Schapery’s  expressions  were  used  to  obtain  the  yarn  properties  from  the  given  fiber  and 
matrix  properties.  Tht^iTEx-\0  results  for  examples  6 and  7 were  compared  with  TEXCAD 
— an  approximate  analytical  method  developed  by  Naik  (1994).  In  both  the  examples  the 
overall  fiber  volume  fraction  was  0.64.  The  unit-cell  was  divided  into  13x13x7  elements  for 
all  three  examples.  The  elastic  constants  for  the  three  textile  composite  examples  are 
presented  in  Table  5.  It  must  be  noted  that  jj,TEx—\Q  will  marginally  under-predict  the 
stiffness  moduli — since  the  yam  cross-section  in  the  numerical  model  is  approximated  as  a 
polygon  inscribed  within  the  actual  cross-sectional  area.  Consequently  the  yarn/fiber 
volume  fraction  in  the  numerical  model  will  always  be  less  than  the  theoretical  volume 
fraction. 

2.2  Stress  Gradient  Effects 


The  methods  explained  in  Section  2. 1 assume  that  the  unit-cells  exist  in  all  the  three 
directions.  This  will  be  true  in  the  case  of  thick  textile  composites.  However  there  are  many 
applications  in  which  thin  composites  are  used.  In  fact  in  order  to  take  advantage  of  the 
properties  of  composites,  the  structures  have  to  be  made  of  thin  plate  like  members  with 
stiffeners  for  load  transfer.  In  such  cases  there  will  be  fewer  unit-cells  in  tbe  thickness 
direction.  Thus  the  free  surface  effects  will  be  predominant.  There  will  be  severe  stress 
gradients  through  the  thickness,  and  they  will  have  an  influence  on  the  apparent  stiffness  and 
strength  of  the  structure  (Marrey  and  Sankar,  1993a;  Sankar  and  Marrey;  1993a). 

The  following  simple  example  will  illustrate  the  stress  gradient  effects  on  stiffness. 
Consider  a layered  medium  consisting  of  alternating  layers  of  materials  of  equal  thickness 
with  Young’s  moduli  E]  and  £2  respectively  (Fig.  12a).  Any  micromechanical  model  would 
predict  that  the  medium  can  be  considered  as  a homogeneous  orthotropic  material  at 
macroscale  and  also  the  effective  Young’s  modulus  in  the  longitudinal  direction  is 
(E]  +E2)/2,  and  there  would  not  be  any  bending-stretching  coupling  in  the  principal  material 
direction.  However,  if  we  consider  a bimaterial  beam  consisting  of  the  same  two  materials 
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(Fig.  12b),  we  will  find  that  there  is  a bending-stretching  coupling,  and  also  the  flexural 
rigidity  cannot  be  predicted  from  the  Young’s  modulus  of  the  homogeneous  orthotropic 
medium  and  the  total  beam  thickness.  The  bimaterial  beam  has  properties  and  behavior 
different  from  the  corresponding  infinite  medium.  This  phenomenon  is  observed  in  the 
transverse  shear  behavior  also  (Sankar  and  Marrey,  1993b).  A similar  behavior  is  also 
expected  in  thin  textile  composites  where  there  are  fewer  unit-cells  in  the  thickness 
directions,  and  the  unit-cells  are  not  subjected  to  a macroscopically  homogeneous  state  of 
deformation  as  assumed  in  Section  2.1. 


One  method  of  overcoming  the  above  difficulties  in  thin  textile  composites  is  to 
model  the  composite  as  a plate/beam,  and  compute  the  structural  stiffness  properties  (e.g., 
[A],  [5]  and  [D]  of  the  plate)  directly  from  the  unit-cell  analysis  instead  of  the  continuum 
stiffness  properties  such  as  Young’s  modulus,  shear  modulus  etc.  In  the  following  sections 
we  illustrate  these  concepts — first  for  a thin  textile  composite  modeled  as  a beam  (Sankar 
and  Marrey,  1993b;  Marrey  and  Sankar,  1993b)  and  then  for  a textile  composite  plate.  The 
purpose  of  the  beam  model  is  to  present  the  issues  involved  in  computing  the  structural 
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stiffness  coefficients  as  opposed  to  continuum  elastic  constants.  Further  the  periodic  EC’s 
are  different  from  those  in  the  continuum  model. 


2.3  Unit-Cell  Analysis  for  Beam  Thermo-Mechanical  Coefficients 


In  this  section,  a procedure  for  finding  the  equivalent  flexural  stiffness  properties  and 
thermal  coefficients  of  a textile  structural  composite  beam  is  described.  The  textile 
composite  beam  is  assumed  to  be  in  the  xz-plane  with  unit-cells  repeating  in  the  x-direction. 
A state  of  plane  strain  parallel  to  the  xz-plane  is  assumed.  On  the  macroscale  it  is  assumed 
that  the  beam  is  homogeneous  and  its  behavior  can  be  characterized  by  the  following  beam 
constitutive  relation  : 
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where  P,  M and  V are  the  axial  force,  bending  moment  and  transverse  shear  force  resultants, 
respectively;  [K\  is  the  symmetric  matrix  of  beam  stiffness  coefficients;  ?<^andyo^are 

the  midplane  axial  strain,  curvature  and  transverse  shear  strain  respectively;  ap,  and  ay, 
respectively,  are  the  beam  thermal  expansion,  thermal  bending  and  thermal  shear 
coefficients.  The  midplane  deformations  are  related  to  the  midplane  axial  displacement  Uo, 
transverse  displacement  w,  and  rotation  xp  as: 
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Actually,  A77,  Kj2,  K22  and  ^55  are  similar  to  the  laminate  stiffness  coefficients  A77,  Bj],  Du 
and  respectively.  There  is  no  equivalence  for  Kj^  and  K23  in  the  laminate  theory, 

because  the  layers  are  assumed  to  be  orthotropic,  and  they  are  rotated  about  the  z-axis  only. 
However,  such  a coupling  between  inplane  deformations  and  transverse  shear  deformation 
may  exist  in  textile  composites  as  the  fibers  are  inclined  to  the  xy-plane  unlike  in  the 
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laminates.  The  beam  constitutive  relation  in  Eqn.  (12)  can  also  be  expressed  in  terms  of 
compliance  coefficients  as: 
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2.3.1  Beam  Steady  State  Loading  Conditions 


(d) 


Figure  13.  Steady  state  loading  for  a beam. 


The  continuum  unit-cell  analysis  assumes  that  all  the  unit-cells  are  subjected  to 
identical  stress  and  strain  fields,  for  a given  state  of  loading.  This  is  true  in  the  case  of 
constant  axial  force  (P)  and  constant  bending  moment  (Af)  in  the  beam  (Figs.  13a  and  13b). 


27 


However,  when  a shear  force  (V)  is  applied  to  the  beam,  the  shear  force  will  give  rise  to 
building  up  of  bending  moment  at  every  cross-section,  such  that  V=-  (dM/dx).  This  bending 
moment  varies  linearly  over  the  length  of  the  beam  violating  the  assumption  of 
homogeneous  deformation.  A state  where  the  unit-cells  are  subjected  to  identical 
deformation  under  a shear  force  can  be  created  by  adding  a couple  periodically  (Fig.  13c), 
or  by  having  shear  tractions  on  the  top  and  bottom  surfaces  to  cancel  the  bending  moment 
continuously  (Fig.  13d).  In  both  cases  traction-free  conditions  are  violated  on  the  top  and 
bottom  surfaces  of  the  beam.  As  will  be  seen  later,  this  situation  creates  difficulties  in 
estimating  the  shear  stiffness  of  the  beam  accurately. 

2.3.2  Unit-Cell  Boundary  Conditions 

The  left  end  of  the  unit  cell,  x=0  (Fig.  14)  are  subject  to  minimum  support  constraints 
to  prevent  rigid  body  translation  and  rotation.  The  top  and  bottom  surfaces  of  the  beam  are 
assumed  to  be  traction  free.  The  edges  x=0  and  x=L  have  identical  nodes  in  the  finite  element 
model  and  periodic  EC’s  are  enforced  between  these  nodes.  Three  linearly  independent 
deformations  were  applied  to  the  unit-cell,  in  order  to  find  the  beam  stiffness  matrix  [/f] 
namely. 

Case  (i)  unit  axial  strain  = 0,y^  = 0) 

Case  (ii)  unit  curvature  accompanied  by  transverse  deflection  such  that  the  transverse 

shear  strain  was  zero  (e^  = 0,x^  = = 0) 

Case  (iii)  unit  transverse  shear  strain  (Eo  = 0,x’^  = 0,y^  = 1) 

The  periodic  displacement  constraints  applied  for  each  deformation  case  are  presented  in 


Table  6. 
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Table  6.  Periodic  displacement  EC’s  for  beam  unit-cell. 


u(L,z)-u(0,z) 

w(L,z)-w(0,z) 

AT 

Case  i.  unit  axial  strain  (so^=l) 

L 

0 

0 

Case  ii.  unit  curvature  (x^=l) 

Lz 

-L^/2 

0 

Case  iii.  unit  shear  strain  (yq^=1  ) 

0 

L 

0 

Consider  the  second  deformation  case,  where  a curvature  is  applied  to  the  unit-cell.  Let ’m’ 
and  ’n’  be  a set  of  corresponding  nodes  on  the  left  and  right  ends  of  the  unit-cell.  The  periodic 
EC’s  imposed  between  the  nodes ’m’  and  ’n’  are  shown  in  Fig.  14.  Applying  a curvature  to 
the  unit-cell  also  induces  transverse  shear  strain  due  to  the  V'^term  in  the  expression  for 
(Eqn.  13).  The  difference  in  displacement  in  the  z-direction  is  therefore  imposed  to  make  the 
macroscopic  transverse  shear  strain  equal  to  zero. 


z 


F"*  = — F" 
^ z ^ z 


Figure  14.  Eoundary  conditions  pertaining  to  the  deformation  to  0,  So^=0, 

Yo^=0. 
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2.3.3  Determination  of  Beam  Stiffness  Coefficients 


The  three  linearly  independent  deformations  described  in  the  preceding  section  are 
applied  to  the  unit-cell.  The  temperature  difference  {AT)  is  set  to  zero  for  all  three 
deformations.  For  each  case,  the  axial  force  P,  the  bending  moment  Me  and  the  shear  force 
resultant  V are  computed  from  the  nodal  forces  at  the  ends  of  the  unit-cell.  As  explained 
earlier,  the  axial  force  and  shear  force  are  constant  through  the  unit-cell,  but  the  bending 
moment  varies  linearly  across  the  unit-cell  length.  Therefore  the  bending  moment  at  the 
center  of  the  unit-cell  is  used  as  the  reference  moment,  which  is  obtained  by  averaging  the 
bending  moment  resultants  at  the  ends  of  the  unit-cell.  From  the  force  resultants,  a 
pseudo-stiffness  matrix  [/:]  can  be  computed  that  relates  the  forces  and  deformations  as 


^11  ^12  ^13 
^21  ^22  ^23 
^31  ^32  ^33 

For  example,  ku,  k2i  and  ksj  will  correspond  to  the  values  of  P,  Me  and  V obtained  for  the 
case  of  unit  axial  extension  of  the  unit-cell  (Case  i).  The  pseudo  stiffness  matrix,  in  general 
will  not  be  symmetric  since  it  does  not  relate  corresponding  forces  and  deformations 
(conjugate  quantities,  product  of  which  yields  an  energy  term) — it  can  be  rather  considered 
as  a matrix  of  influence  coefficients.  The  inverse  of  [^] , denoted  by  [i]  has  some  significance. 
The  matrix  [s]  is  defined  as 
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We  know  that  V=0  is  a steady  state  loading  condition  such  that  Me=M,  where  M is  the 
constant  bending  moment  along  the  unit-cell.  Noting  that  only  the  last  column  of  [s] 
multiplies  with  V,  and  comparing  Eqns.  (14)  and  (16),  we  can  conclude  that  the  first  two 
columns  of  [5]  and  [5]  must  be  identical  to  each  other.  Since  [5]  is  symmetric  (Si3=S3], 
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^23=^32),  we  can  determine  all  but  S33  of  the  beam  compliance  matrix  [5].  An  alternative 
approach  to  predict  the  shear  compliance  S33  is  explained  in  the  next  section. 

2.3.4  Determination  of  Beam  Shear  Stiffness 

The  difficulty  in  estimating  S3 3 (or  K3 3)  is  associated  with  the  inability  to  create  a 
state  of  deformation  such  that  only  V is  present.  The  shear  modulus  of  the  textile  composite 
Gxz  can  be  computed  by  assuming  that  the  unit-cells  span  the  material  in  the  z-direction  also. 
In  that  case  periodic  EC’s  are  imposed  between  the  top  and  bottom  surfaces  of  the  unit-cell, 
and  the  unit-cell  is  deformed  transversely.  Then,  there  will  be  shear  tractions  on  the  top  and 
bottom  surfaces  of  the  beam — in  fact  this  situation  would  correspond  to  Fig.  1 3(d).  One  may 
surmise  that  a shear  correction  factor  could  be  found  such  that  K33=x^Gxzh.  But  a simple 

bimaterial  beam  example  will  show  that  the  shear  stiffness  can  be  grossly  underestimated. 

Consider  a bimaterial  beam  with  layers  of  equal  thickness  (h/2).  By  performing  a 
continuum  unit-cell  analysis  the  shear  modulus  of  the  material  is  found  to  be  equal  to  2Gj  G2 
/ (Gy-f-  G2).  Assuming  G;/G2=  10,  Gxz=  0.182  Gj.  However  the  actual  shear  stiffness  K33 
is  equal  to  Gj  -1-  G2)  h/2,  which  means  that  the  apparent  shear  stiffness  Gxz=  (Gj+  G2)/2 
= 0.55G] . This  is  about  three  times  the  previous  estimate.  This  discrepancy  is  due  to  differing 
assumptions  regarding  the  constancy  of  shear  stress  or  shear  strain.  The  continuum  unit-cell 
analysis  imposes  constant  shear  stress  in  the  two  materials,  and  hence  the  shear  compliance 
of  the  composite  is  the  average  of  the  compliances  of  the  constituent  materials.  This  is  true 
when  there  are  a large  number  of  unit-cells  in  the  z-direction.  In  a bimaterial  beam,  however, 
the  shear  strain  is  almost  constant  in  the  two  layers,  and  hence  the  shear  stiffness  is  the 
average  of  the  shear  stiffness  of  the  individual  layers — which  is  consistent  with  the  method 
of  computing  A55  in  the  lamination  theory.  This  illustrates  the  need  for  special  procedures 
for  predicting  the  shear  stiffness  of  thin  textile  composite  beams. 

To  overcome  this  problem  we  use  an  energy  approach  to  compute  the  shear  stiffness 
of  a thin  textile  composite.  We  model  a beam  of  length  2L  consisting  of  two  unit-cells,  with 
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plane  strain  finite  elements.  The  beam  is  subject  to  boundary  conditions  corresponding  to 
pure  shear  strain  (third  boundary  condition  in  Table  6).  The  top  and  bottom  surfaces  of  the 
beam  are  traction  free.  The  shear  strain  energy  over  a length  L in  the  middle  of  the  beam,  U^, 
is  computed  from  the  finite  element  results  as 


where  and  the  shear  stress  and  shear  strain  at  the  center  of  the  element  and 

A,-  is  the  area  of  the  element.  The  summation  is  performed  over  all  elements  located  in  a 

length  L in  the  middle  of  the  beam.  Next,  the  shear  strain  energy  over  the  same  length  L is 
computed  using  the  beam  formula: 

Us  = + >^33^^  (18) 

In  the  above  equation  P,  Me  and  V can  be  obtained  from  the  finite  element  results.  The 
coefficients  5/j  and  S23  have  already  been  estimated.  Then  by  equating  the  shear  strain 
energy  quantities  in  Eqns.  (17)  and  (18),  the  only  unknown,  that  is  S33,  can  be  evaluated.  The 
choice  of  two  unit-cells  to  perform  the  above  analysis  deserves  an  explanation.  When  this 
was  tried  with  one  unit-cell  for  the  cases  of  isotropic  and  bimaterial  beam,  the  results  were 
not  good.  It  was  mentioned  earlier,  that  the  application  of  a shear  force  would  result  in  a 
couple  at  the  unit-cell  ends  to  continuously  cancel  the  effect  of  the  bending  moment  created 
along  the  beam  length.  This  couple  is  manifested  in  the  finite  element  results  as  concentrated 
forces  at  the  four  corners  of  the  unit-cell,  thus  creating  severe  stress  concentrations.  When 
two  unit-cells  are  used  in  the  model,  the  stress  concentrations  remain  in  the  corners  of  the 
beam,  but  their  effects  diminish  in  the  middle  portion  of  the  beam.  As  will  be  seen  in  the  later, 
the  two  unit-cell  method  gave  very  good  K33  for  both  isotropic  and  bimaterial  beams. 

2.3.5  Determination  of  Beam  CTE’s 

The  procedure  for  determining  the  beam  thermal  expansion  coefficients  (Marrey  and 
Sankar,  1993b)  is  as  follows.  The  beam  unit-cell  is  subject  to  a uniform  temperature 
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difference  given  hy  AT  = Tq.  The  deformations  in  the  beam  are  restrained  by  setting 
£o’^=‘>i^=yo^=0.  Then  the  beam  constitutive  relation,  Eqn.  (12),  will  reduce  to 
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The  axial  force  P,  the  bending  moment  resultant  at  the  center  of  the  unit  cell  Me  and  the  shear 
force  V are  computed  from  the  nodal  forces  at  the  ends  of  the  unit-cell.  Then  the  beam  CTE’s 
can  be  estimated  from  the  expression. 
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2.3.6  Results  and  Discussion 


The  procedures  described  above  were  implemented  for  the  following  cases: 

(a)  an  isotropic  beam 

(b)  a bimaterial  beam  with  isotropic  layers  of  equal  thickness 

(c)  a plain  weave  textile  composite  beam  where  the  yam  is  assumed  to  be  transversely 
isotropic  and  the  matrix  is  isotropic. 


Table  7.  Constituent  material  properties  for  beam  examples. 


isotropic  beam 

E=10  GPa,  V = 0.30,  a = lOxlO-^  PC 

bimaterial  beam 

Ej  = 70  GPa,  V]  = 0.33,  aj  = 23x10“^  PC 
E2  = 3.5  GPa,  V2  = 0.35,  «2  = 60x10-^  PC 

plain-weave 
textile  beam 

yarn: 

E;  = 159  GPa,  £2  = 10.9  GPa,  G/2  = 6.4  GPa, 
vj2=  0.38,V2j  = 0.38,  = 0.045x10-6  PC,  a2  = 20.2x1 0"6  PC 

where  the  yam  direction  is  parallel  to  the  i-axis  and  25-plane 
is  the  plane  of  isotropy. 

matrix: 

Err,  = 3.5  GPa,  Vm  = 0.35,  am  = 60x10-6  /oq 
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The  properties  of  the  constituent  materials  for  all  the  cases  are  listed  in  Table  7.  The 
dimensions  of  the  unit-cell  and  the  yarn  architecture  for  the  textile  beam  were  taken  from 
Yoshino  and  Ohtsuka  (1982).  The  same  unit-cell  dimensions  (length  of  3.6  mm  and  height 
1.8  mm)  were  also  used  for  the  isotropic  and  bimaterial  cases.  The  unit-cell  of  the  beam  was 
discretized  using  eight-node  isoparametric  plane  strain  finite  elements.  The  finite  element 
mesh  for  the  isotropic  unit-cell  and  the  plain  weave  unit-cell  were  identical  except  that 
different  material  properties  were  used. 


Table  8.  Comparison  of  beam  stiffness  coefficients  and  CTE’s  (SI  units). 


Kii 

Ki2 

K22 

K33 

ap  PC 

aMPC  II 

isotropic 

beam 

unit-cell 

analysis 

19.78x106 

0 

5.35 

5.96x106 

10x10-6 

0 

beam 

theory 

19.78x106 

0 

5.34 

5.77x106 

lOxKH 

0 

bimaterial 

beam 

unit-cell 

analysis 

74.29x106 

30.20x103 

20.06 

8.47x106 

30.73x10^ 

-14.62x10-3 

beam 

theory 

74.29x106 

30.20x103 

20.06 

8.62x106 

30.74x10^ 

-14.63x10-3 

textile 

beam 

unit-cell 

analysis 

27.76x106 

0 

5.41 

9.21x106 

12.66x10-^ 

-24.12x10-^ 

mosaic 

model 

71.48xl(H 

0 

8.13 

8.14x106 

4.39x10-6 

0 

Note.-  Kj3,  K23  and  ay  are  zero  for  all  cases 


The  deformed  unit-cells  under  various  independent  loading  conditions  are  shown  in 
Figs.  15  to  17.  The  stiffness  and  thermal  coefficients  for  the  three  beams  are  shown  in  Table 
8.  The  results  for  the  isotropic  and  bimaterial  beams  were  compared  to  exact  beam  theory 
solutions.  Exact  shear  correction  factors — 0.833  for  the  isotropic  beam  and  0.555  for 
bimaterial  beam  (Whitney,  1973)  were  used  in  the  beam  theory  solution  to  compute  the  shear 
stiffness  given  in  Table  8.  It  can  be  seen  from  Table  8 that  the  beam  unit-cell  analysis  is  able 
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to  predict  the  axial  and  bending  stiffness  coefficients  {Kjj  and  K22)  very  accurately.  As 
expected  the  shear  stiffness  {K33  or  A55)  predictions  have  errors,  but  they  are  very  minimal. 
It  can  be  noticed  that  the  comers  of  the  unit-cell  are  severely  deformed  (Figs.  15c,  16c  and 
17d),  when  the  unit-cell  is  subject  to  constant  shear  strain  leaving  the  top  and  bottom  surfaces 
traction  free.  However  when  shear  tractions  are  allowed  on  the  top  and  bottom  surfaces  of 
the  unit-cell,  the  distortions  at  the  corners  disappear  (Figs.  15d,  16d  and  17e).  Then  what  is 
obtained  is  the  shear  modulus  Gxz  and  not  the  beam  shear  stiffness.  The  shear  modulus  of 
the  plain  weave  beam  was  found  to  be  3.07  GPa.  This  would  yield  the  apparent  shear  stiffness 
asGxzh=5.53x  10^Nm~' ; whereas  the  actual  shear  stiffness  is  9.21  x 10^Nm~*  (K33  in  Table 
8).  The  Young’s  modulus  of  the  textile  beam  Ex  may  be  extracted  from  Kj],  as  Kji/h,  which 
would  yield  Ex  = 15.42  GPa.  If  this  value  of  Ex  were  used  to  predict  the  flexural  stiffness  of 
a homogeneous  beam  as  Djj  = Exh^/12,  we  would  obtain  D77  as  7.50  Nm — whereas  the 
actual  flexural  stiffness  is  equal  to  5.41  Nm.  The  same  idea  holds  for  the  beam  thermal 
coefficients  also.  The  beam  CTE’s  ap,  and  ay  cannot  be  predicted  from  the 
corresponding  continuum  CTE’s.  Table  9 shows  the  discrepancy,  for  the  plain  weave 
example,  between  the  beam  CTE’s  obtained  directly  and  the  beam  CTE’s  obtained  from  the 
corresponding  continuum  CTE’s.  It  may  be  noted  that  the  continuum  model  would  always 
predict  the  thermal  expansion  coefficient  ap  as  ax,  and  the  thermal  bending  coefficient 
as  zero.  This  underscores  the  importance  of  the  present  analysis  for  predicting  the  beam 
stiffness  properties  for  a thin  textile  composite  directly. 


Table  9.  Comparison  of  textile  CTE’s. 


CTE’s  from 

CTE’s  from 

% error 

beam  model 

continuum  model 

ap  X 10-6  /t>C 

12.66 

11.30 

-10.46 

aM  X 10-6  PC/m 

-24.12 

0 

00 

ay/^C 

0 

0 

0 
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(a) 


(b) 


(c) 


(d) 


Figure  15.  Deformed  unit-cell  of  the  isotropic  beam,  (not  to  scale) 

(a)  unit  extensional  strain;  (b)  unit  curvature;  (c)  unit  shear  strain,  top  and 
bottom  surfaces  are  traction  free;  (d)  unit  shear  strain,  tractions  allowed  on 
top  and  bottom  surfaces. 
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Figure  16.  Deformed  unit-cell  of  the  bimaterial  beam,  (not  to  scale) 

(a)  unit  extensional  strain;  (b)  unit  curvature;  (c)  unit  shear  strain,  top  and 
bottom  surfaces  are  traction  free;  (d)  unit  shear  strain,  tractions  allowed  on 
top  and  bottom  surfaces. 


37 


Figure  17.  Textile  beam,  (not  to  scale) 

(a)  undeformed  unit-cell;  and  deformation  under:  (b)  unit  extensional  strain; 
(c)  unit  curvature;  (d)  unit  shear  strain,  top  and  bottom  surfaces  are  traction 
free;  (e)  unit  shear  strain,  tractions  allowed  on  top  and  bottom  surfaces. 
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The  textile  beam  stiffness  coefficients  were  also  estimated  using  a procedure  similar 
to  the  mosaic  model  (Ishikawa  and  Chou,  1982a).  They  are  compared  with  the  coefficients 
obtained  from  the  unit-cell  analysis  (Table  8).  The  idealization  made  in  the  mosaic  model 
is  shown  in  Figure  18.  The  unit  cell  is  divided  into  five  segments,  with  each  segment  modeled 
as  a cross-ply  laminate  with  a stacking  sequence  which  best  represented  the  yarn  architecture 
within  that  segment.  The  stiffness  matrix  of  each  segment  was  computed  using  laminate 
analysis.  Then  the  compliance  of  the  textile  beam  was  computed  as  the  length-weighted 
average  of  the  compliance  of  the  five  segments. 

5 

[S]  = (^)  X (21) 

^ k=\ 

From  Table  8,  it  can  be  seen  that  the  mosaic  model  predicts  K33  reasonably  well.  The  reason 
for  the  lack  of  agreement  in  Kjj  and  K22  can  be  attributed  to  the  fact  that  a major  portion  of 
the  yarn  is  modeled  as  a 0°  laminate  in  the  mosaic  model,  which  tends  to  over-predict  the 
axial  and  flexural  stiffnesses. 


Figure  18.  Mosaic  model  of  the  textile  unit-cell. 
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2.4  Plate  Thermo-Mechanical  Coefficients 


The  previous  section  described  the  concepts  involved  in  modeling  a thin  textile 
composite  structure  as  as  homogeneous  beam  to  predict  the  beam  stiffness  coefficients  and 
CTE’s.  The  procedure  is  extended  in  this  section,  utilizing  three-dimensional  finite  element 
analysis,  to  model  the  thin  composite  as  a plate  (Sankar  and  Marrey,  1992)  to  determine  the 
corresponding  thermo-elastic  coefficients. 

The  plate  is  assumed  to  be  in  the  xy-plane  with  unit-cells  repeating  in  the  x-  and 
y-directions.  The  lengths  of  the  unit-cell  in  the  x-  and  y-directions  are  assumed  to  be  a and 
b respectively  and  the  unit-cell  thickness  as  h.  On  the  macroscale  the  plate  is  assumed  to  be 
homogeneous  and  the  plate  behavior  is  characterized  by  the  plate  constitutive  relation: 
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My 
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xy 
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where  £,o^,  and  are  the  midplane  axial  strain,  shear  strain  and  curvature;  a^P  and 
PiP  are  the  plate  thermal  expansion  and  bending  coefficients;  A/,-  and  Mi  are  the  axial  force 
and  bending  moment  resultants  respectively  in  the  homogeneous  plate.  The  plate  stiffness 
matrix  comprises  of  the  [A],  [B]  and  [D]  sub-matrices,  which  are  the  plate  extensional 
stiffness,  coupling  stiffness  and  bending  stiffness  matrices  respectively.  The  plate  stiffness 
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matrix  can  be  expanded  as: 
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fi^l2 

fi^22 

Die 

^16 

^26 
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The  midplane  strains  and  curvatures  are  related  to  the  midplane  displacements  and  rotations 
as: 
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(24) 

(25) 


2.4.1  Unit-Cell  Boundary  Conditions 

The  plate  thermo-mechanical  properties  are  obtained  by  modeling  the  unit-cell  with 
eight-node  brick  elements  and  subjecting  the  unit-cell  to  six  linearly  independent 
deformations.  The  six  deformations  are  given  by:  (1)  unit  maintaining  the  rest  of  the 
macroscopic  strains  and  curvatures  as  zero;  (2)  unit  such  that  remaining  strains  and 
curvatures  are  zero;  and  similarly  (3)  unit  y^;  (4)  unit  (5)  unit  (6)  unit  x^.  In  the 
last  three  cases  (non-zero  curvatures)  the  deformation  was  accompanied  with  a transverse 
deflection  such  that  the  transverse  shear  strain  was  zero  (Table  10). 
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Table  10.  Periodic  displacement  EC’s  imposed  on  the  lateral  faces  of  the  plate  unit-cell. 


u(a,y)~ 

u(0,y) 

v(a,yh 

v(O.y) 

wfa.yj- 

w(0,y) 

u(x,b)- 

u(x,0) 

v(x.b)~ 

v(x,0) 

w(x,b)- 

w(x,0) 

AT 

1. 

Sx0^=l 

a 

0 

0 

0 

0 

0 

0 

2. 

Cy0^=l 

0 

0 

0 

0 

b 

0 

0 

3. 

II 

0 

a/2 

0 

b/2 

0 

0 

0 

4. 

II 

az 

0 

-a^/2 

0 

0 

0 

0 

5. 

0 

0 

0 

0 

bz 

-b^/2 

0 

6. 

Xxy^=l 

0 

azJ2 

-ay/2 

bz/2 

0 

-bx/2 

0 

7. 

AT^=1 

0 

0 

0 

0 

0 

0 

1 

b 

^v=0  1 

■ H=0,  v=0 

U=0,  w=0 
X 

Figure  19.  Eoundary  conditions  on  plate  unit-cell  to  restrain  rigid  body 
translation  and  rotation. 

The  unit-cell  is  subject  to  minimum  support  constraints  to  prevent  rigid  body  rotation 
and  translation  (Fig.  19).  The  top  and  bottom  surfaces  of  the  plate  are  assumed  to  be  free  of 
tractions.  The  faces  x=0  and  x=a  have  identical  nodes  in  the  finite  element  model,  and  so 
do  the  pair  of  faces  y=0  and  y=b.  The  identical  nodes  on  opposite  faces  of  the  unit-cell  are 
constrained  to  enforce  the  periodic  EC’s.  The  traction  boundary  conditions  on  the  lateral 
faces  of  the  unit-cell  are: 
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Fi{a,y,z)  = - FiiQ,y,z) 

F^{x,b,z)  = - Flx,Q,z) 


(26) 


The  periodic  displacement  EC’s  enforced  for  each  mode  of  deformation  are  presented  in 
Table  10. 


2.4.2  Determination  of  Plate  Stiffness  Coefficients  and  CTE’s 


The  six  linearly  independent  deformations  are  applied  to  the  unit-cell  such  that  only 
one  of  the  macroscopic  strains  or  curvatures  is  non-zero  (first  six  cases  in  Table  10).  The 
temperature  difference  is  set  to  zero  for  all  six  cases.  It  must  be  noted  that  the  applied 
deformations  must  ensure  that  the  transverse  shear  strains,  yxz’^  and  are  zero  where 


(27) 


The  force  and  moment  resultants  can  be  obtained  by  one  of  the  following  two  methods.  In 
the  first  method,  the  resultants  are  computed  by  averaging  the  nodal  forces  on  each  face  of 
the  unit-cell.  For  example,  on  the  face  x=a  the  force  and  moment  resultants  are  computed 
using  the  relations: 


Nx  = 


Nxy  = 


Mx  = 


(28) 


(\  ” 

i)  Ff(a,y,z) 

where  Fx^^^  and  Fy(^^  are  the  nodal  forces  in  the  x and  y directions  at  the  node  and  is 
the  total  number  of  nodes  on  the  face.  The  force  and  moment  resultants  can  also  be  computed 
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by  averaging  the  microstresses  over  the  unit-cell  volume.  Then  the  resultants  on  the  face  x=a 
are  obtained  as 


Nr  = 


±( 

ab  J 


oAx,y,z)dV  Nxy  = 
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ab  J 


Xxy{x,y,z)dV 
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K V 

Substituting  the  values  of  the  deformation  and  the  force  resultants  in  the  plate  constitutive 
relation,  Eqn.  (22),  the  stiffness  coefficients  in  the  column  corresponding  to  the  non-zero 
deformation  can  be  computed.  This  procedure  is  repeated  for  other  deformation  components 
to  obtain  all  the  stiffness  coefficients. 

To  predict  the  CTE’s,  the  plate  unit-cell  is  subject  to  a uniform  temperature 
difference,  given  by  AT  = Tq.  In  the  finite  element  model,  periodic  displacement  EC’s  are 
applied  such  that  all  six  components  of  the  deformation  are  zero  (seventh  case  in  Table  10). 
The  averaged  force  and  moment  resultants  are  computed  using  one  of  the  procedures 
described  above.  The  thermal  expansion  coefficients  aP,  and  thermal  bending  coefficients 
PP  are  then  obtained  from  the  relation: 
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2.4.3  Results  for  Plate  Stiffness  Coefficients 


The  plate  [A],  [5],  [Z)]  matrices  and  CTE’s  were  found  for  the  seven  examples  listed 
in  Table  2 by  implementing  the  finite  element  code  ^TE;f-10.  As  a preliminary  check,  the 
code  was  executed  for  an  isotropic  plate  and  compared  with  the  plate  properties  using 

Fh 

lamination  theory  (for  one-ply).  For  example,  A77  was  calculated  as  — ■ and  D77  as 

(1  - v2) 


Eh^ 

12(1  - v2) 


•.  Then  the  properties  were  computed  for  a bimaterial  plate  using  the  finite 
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element  code.  The  plate  properties  for  the  isotropic  and  bimaterial  cases  are  presented  in 
Tables  1 1 and  12  respectively.  The  bimaterial  plate  properties  were  also  computed  using  the 
lamination  theory  for  two  plies,  and  from  the  homogenous  3-D  elastic  constants  computed 
in  the  previous  section  (Section  2.1).  For  example  the  coefficient  Djj  is  obtained  from  the 


3-D  elastic  constants  as  Djj  = 


12(1  - vp 


. The  finite  element  results  for  the  bimaterial 


case  were  exact,  i.e.,  identical  to  the  results  obtained  with  the  two-ply  lamination  theory.  The 
[D]  matrix  computed  from  the  bimaterial  3-D  constants  was  found  to  be  in  good  agreement 
with  the  two-ply  lamination  theory  only  because  both  the  layers  were  equal  in  thickness.  This 
is  a special  case,  and  in  general,  the  [D]  matrix  obtained  from  the  3-D  elastic  constants  will 
be  different  from  the  two-ply  lamination  theory  results. 

The  plate  properties  for  the  unidirectional  composite  examples  are  presented  in  Table 
13  and  for  the  textile  examples  in  Table  14.  In  all  the  examples  it  was  found  that  the  plate 
properties,  especially  , {aP}  and  {^p]  could  not  be  predicted  from  the  corresponding 

3-D  elastic  constants. 


Table  11.  Non-zero  [A],  [B]  and  [D]  coefficients  for  example  1 (isotropic  plate)  using 
finite  elements 


An 

xlO^ 

An 

xlO^ 

A22 

xlO^ 

A^6 

xlO^ 

Dll 

xlO-3 

Di2 

xlO-3 

D22 

xlO-3 

D66 

xlO-3 

GxP,  GyP 

xitr^/oc 

^TEx-\0 

(FEA) 

2.810 

0.843 

2.810 

0.983 

15.320 

4.606 

15.320 

5.358 

10 

lamination 

theory 

2.810 

0.843 

2.810 

0.983 

15.310 

4.593 

15.310 

5.358 

10 

Note:  [A],  [S]  and  [D]  coefficients  in  SI  units 
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Table  12.  Non-zero  [A],  [j5]  and  [D]  coefficients  for  example  2 (bimaterial  plate)  using 
finite  elements 


^Ih  ^22 
xl06 

An 

x1Q6 

A66 

xl06 

Blh  ^22 
xlO^ 

Bn 

xl03 

^66 

xl0~^ 

fiTEx-m 

(FEA) 

9.832 

2.043 

3.895 

- 0.563 

-0.108 

- 0.228 

lamination  theory 
for  two  plies 

9.832 

2.043 

3.895 

-0.563 

-0.108 

- 0.228 

lamination  theory 
using  3-D  elastic 
constants 

9.844 

2.048 

3.899 

0 

0 

0 

E>ii,  D22 

Dn 

B>66 

ttjcP,  ayP 

xl0~^ 

xl0“^ 

xl0“^ 

xlO-6/oC 

PCIm 

fiTEx~\0 

(FEA) 

53.590 

11.149 

21.220 

17.800 

0.170 

lamination  theory 
for  two  plies 

53.573 

11.131 

21.220 

17.814 

0.170 

lamination  theory 
using  3-D  elastic 
constants 

53.762 

11.183 

21.293 

8.190 

0 

Note:  [A],  [5]  and  [D]  coefficients  in  SI  units 
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Table  13.  Non-zero  [A],  [5]  and  [D]  coefficients  for  examples  3 and  4 (unidirectional 
composite)  using  finite  elements 


An 

xlO^ 

An 

xl06 

A22 

xl06 

A66 

Xl06 

Example  3 

^TEx-lO 

(FEA) 

0.690 

0.149 

0.496 

0.177 

Halpin-Tsai  Eqns. 
and  lamination 
theory 

0.673 

0.109 

0.363 

0.113 

Example  4 

juTEx-lO 

(FEA) 

0.452 

0.062 

0.285 

0.114 

Halpin-Tsai  Eqns. 
and  lamination 
theory 

0.444 

0.039 

0.151 

0.045 

Dn 

xlO-6 

Dn 

xlO-6 

D22 

xlO-^ 

D66 

xlO-^ 

axPx 

10-^ 

ttyPX 

10-6 

Example  3 

fxTEx~\0 

(FEA) 

3.589 

0.596 

1.980 

0.947 

15.489 

26.184 

Halpin-Tsai 
Eqns.  and 
lamination 
theory 

5.606 

0.908 

3.026 

0.939 

15.625 

55.112 

Example  4 

^lTEx-\0 

(FEA) 

2.256 

0.224 

0.873 

0.568 

131% 

13.188 

Halpin-Tsai 
Eqns.  and 
lamination 
theory 

3.702 

0.328 

1.262 

0.371 

6.774 

34.239 

Note:  [A],  [B]  and  [D]  coefficients  in  SI  units 
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Table  14.  Non-zero  [A],  [5]  and  [D]  coefficients  for  examples  5,  6 and  7 using  finite 
elements 


^11 . ^22 
x1Q6 

An 

xl06 

A66 

xl06 

Bii  xl03 

Example  5 

fiTEx-\0 

(FEA) 

2.681 

0.565 

0.489 

0 

lamination  theory 
using  3-D  constants 

2.783 

0.503 

0.490 

0 

Example  6 

fxTEx-\0 

(FEA) 

12.090 

3.470 

1.223 

0 

lamination  theory 
using  3-D  constants 

13.938 

1.787 

1.208 

0 

Example  7 

fxTEx-\0 

(FEA) 

14.683 

1.351 

1.210 

0.495* 

lamination  theory 
using  3-D  constants 

16.531 

0.770 

1.239 

0 

Dji,  D22 
xlO-3 

Dn 

xlO-3 

D66 

xlO-3 

axP,  ttyP 

xlO-6/oc 

PxP 

PCha 

Example  5 

fiTEx-\0 

(FEA) 

5.687 

1.518 

1.577 

27.465 

0 

lamination 
theory  using 
3-D  constants 

12.054 

2.177 

2.124 

28.363 

0 

Example  6 

fA.TEx-\0 

(FEA) 

41.695 

0.373 

5.879 

1.480 

0 

lamination 
theory  using 
3-D  constants 

75.942 

9.734 

6.582 

1.556 

0 

Example  7 

^lTEx-lQ 

(FEA) 

90.072 

1.123 

6.149 

2.910 

-0.037* 

lamination 
theory  using 
3-D  constants 

87.283 

4.195 

6.753 

1.550 

0 

* In  example  7,  822  = -Bji  and  fiyP  = -(ixP 


Note;  [A],  [5]  and  [D]  coefficients  in  SI  units 


CHAPTER  3 

ANALYTICAL  MODELS  FOR  THERMO-ELASTIC  CONSTANTS 


The  complex  yarn  architectures  in  a textile  composite  make  numerical  modeling  of 
the  unit-cell  extremely  difficult.  Besides,  the  computational  memory  and  run-time 
requirements  for  a detailed  finite  element  analysis  are  enormous.  As  explained  in  the 
introductory  chapter,  there  are  several  parameters  that  can  be  changed  to  alter  the  effective 
composite  properties.  These  parameters  may  be  the  fiber  material  in  the  yarn,  fiber  volume 
fraction  in  the  yarn  (also  called  yarn  packing  density),  overall  fiber  volume  fraction,  preform 
architecture  or  the  matrix  material  properties.  This  emphasizes  the  need  for  simple  analysis 
procedures  to  predict  the  trend  in  variation  of  composite  properties  when  one  of  the 
parameters  is  changed.  These  procedures  will  be  of  use  to  a designer  in  determining  the 
optimum  parameters  for  a certain  application. 

Analytical  methods  are  approximate  because  they  assume  certain  forms  for  the  state 
of  stress  and  strain  in  the  unit-cell.  Averaging  the  stiffnesses  or  compliances  of  the  matrix 
and  the  inclusion  has  long  been  used  to  estimate  the  bounds  of  effective  elastic  properties 
of  the  composite.  Essentially  the  stiffness  averaging  assumes  a state  of  uniform  strain  in  the 
composite  (isostrain),  and  compliance  averaging  assumes  a state  of  uniform  stress  (isostress) 
in  the  matrix  and  inclusion.  In  fact  the  rule  of  mixtures  expressions  for  estimating  the 
effective  properties  of  a unidirectional  composite  is  based  on  such  averaging  sehemes.  Naik 
(1994)  proposed  an  analytical  method  (TEXCAD)  in  which  the  yarns  are  discretized  into 
segments.  Knowing  the  direction  of  the  yam  in  each  segment,  the  segment  stiffness  are 
computed  using  appropriate  transformations.  Then  assuming  a state  of  isostrain,  the 
composite  stiffness  is  obtained  by  volume-averaging  the  yarn-segment  stiffness  and  matrix 
stiffness  in  the  unit-cell.  This  method  seems  to  work,  when  there  is  multi-directional 
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reinforcement  in  the  composite.  However  the  method  fails  for  composites  with  preferential 
yarn  reinforcement.  For  example,  the  transverse  modulus  and  the  inplane  shear  modulus  for 
a unidirectional  composite  cannot  be  estimated  using  the  above  method.  Another  analytical 
model  called  the  mosaic  model,  was  proposed  by  Ishikawa  and  Chou  (1982a).  As  discussed 
earlier,  in  the  mosaic  model  the  yarn  architecture  is  simplified  to  that  of  a cross-ply  laminate 
and  the  lamination  theory  is  used  to  predict  the  composite  properties. 

The  state  of  stress/strain  in  a textile  composite,  subjected  to  a uniform  macrostress, 
is  much  more  complex  than  that  assumed  in  the  above  mentioned  methods.  We  propose  a 
scheme  of  selective  averaging — called  the  selective  averaging  method  (SAM) — in  which 
both  stiffness  and  compliance  coefficients  can  be  averaged  selectively  depending  on  a more 
realistic  assumption  of  either  isostress  or  isostrain. 

3.1  Selective  Averaging  Method  IS  AMI  for  Continuum  Properties 

In  this  section,  we  describe  SAM  for  estimating  the  effective  elastic  constants  and 
CTE’s  for  a textile  composite.  Consider  a rectangular  hexahedron  of  dimensions  axbxc 
as  the  unit-cell.  The  unit-cell  is  discretized  into  slices  on  the  mesoscale,  and  elements  on  the 
microscale  as  shown  in  Fig.  20.  To  distinguish  between  the  macrolevel,  mesolevel  and 
microlevel  properties  in  this  section,  an  over-tilde  is  used  to  denote  the  mesolevel  properties, 
and  a superscript  ”M”  is  used  to  denote  the  macrolevel  properties.  For  example,  [C^],  [Q 
and  [C]  will  represent  the  macrolevel,  mesolevel  and  microlevel  stiffnesses  respectively. 

The  objective  here  is  to  determine  the  coefficients  of  the  effective  stiffness  matrix 
[C^]  as  defined  in  Eqn.  (2).  To  find  the  first  column  of  the  effective  stiffness  matrix,  the 
unit-cell  is  divided  into  slices  (mesolevel)  of  thickness  dx  parallel  to  the  yz-plane  (Fig.  20a). 
Each  slice  is  further  sub-divided  into  elements  (microlevel)  as  shown  in  Figs.  20(b)  and 
20(c).  The  unit-cell  is  subjected  to  a deformation  such  that  all  macrostrains  except 
equal  to  zero  and  Exx^=1-  It  is  assumed  that  the  mesolevel  and  microlevel  strains. 
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corresponding  to  the  zero  macrostrains,  are  negligible.  In  other  words, 

ef  = = £.  = 0 i ^ \ (32) 


Figure  20.  Hierarchy  of  discretization  for  a unit-cell  to  determine  first  and  fifth 
columns  of  effective  stiffness  matrix. 

(a)  unit-cell;  (b)  slice;  (c)  element. 


Assuming  a state  of  isostrain  within  the  slice  {£xxi.x,y,  z)  = e{x) ) the  average  stiffness  of 
a slice  can  be  obtained  as: 


c b 

CjiW  = ^ I I Cy^{x,y,z)  dy  dz  (33) 

z=0y=0 

where  C]](x,y,z)  is  the  element  stiffness  coefficient  referred  to  the  unit-cell  coordinates.  The 
stiffnesses  of  the  slices  are  averaged  on  the  macrolevel  based  on  the  isostress  assumption, 
i.e.,  dxx(x)  = <7^ . Then  the  first  column  of  the  effective  stiffness  matrix  can  be  computed 
using  the  following  two  relations: 
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1 

a 


jc  = 0 


Ci,(x) 


■dx 


(34) 


c b a 

CM 

— C^^{x,y,z)  dx  dy  dz  (i  = 2,...,  6)  (35) 

C^^(x) 

z=0y=0x=0 

A similar  procedure  can  be  implemented  to  determine  the  second  and  third  columns  of  the 
macroscale  stiffness  matrix  [C^]. 

A slightly  different  averaging  scheme  is  used  when  the  unit-cell  is  subjected  to  shear 
strains  on  the  macrolevel.  Consider  the  case  where  the  unit-cell  is  subjected  to  unit  at 
macroscale.  The  unit-cell  is  again  discretized  into  sliees  and  elements  as  shown  in  Fig.  20. 
It  is  assumed  that  all  the  other  components  of  strain  at  the  macrolevel,  mesolevel  and 
microlevels  are  zero.  This  can  be  expressed  as: 


pM  _ 


= 


= 0 


4 


(36) 


where  £4  = We  also  assume  that  the  shear  stress  is  constant  in  a slice  such  that 
Tyz(x,y,  z)  = iyzix).  The  shear  compliance  of  a slice  can  then  be  obtained  by  averaging  the 
shear  compliances  of  all  the  elements  in  the  slice  as: 


1 

C44(x) 


J_ 

be 


c b 


n 

z = 0>'  = 0 


1 

C44(x,y,z) 


dy  dz 


(37) 


The  fourth  column  of  the  stiffness  matrix  Q^'^is  obtained  under  the  assumption  that  all  the 
slices  are  under  a state  of  constant  shear  strain: 


/-'M  

'-'(4  ~ 


1 

abc 


c b a 


III 

z=0y=0x=0 


C^jx) 

c^(x,y,z) 


C,4(x,y,  z) 


dx  dy  dz 


a = 1,...,6)  (38) 


A similar  proeedure  is  used  to  determine  the  fifth  and  sixth  eolumns  of  [C^]. 
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To  determine  the  macroscale  CTE’s,  a uniform  temperature  difference  {AT)  is 
applied  throughout  the  unit-cell.  The  unit-cell  is  constrained  from  expanding  such  that  all 
the  strain  components  on  the  macrolevel  are  zero.  A state  of  isostrain  is  assumed  in  the 
unit-cell,  implying  that  the  mechanical  strain  components  on  the  mesolevel  and  microlevel 
are  also  zero.  This  can  be  expressed  as: 

= e,  = £,•  = 0 i = 1,...,6  (39) 

Then  the  thermal  constitutive  relations  on  the  macrolevel  and  microlevel  will  reduce  to: 


[a^]  = - [C^]  {a^]AT 
{a}  = - [C]{a}AT 


(40) 


The  macrostresses  may  be  computed  by  volume-averaging  the  corresponding  microstress 
component  as  shown  below: 


c b a 


±j \ j 


{a}  dx  dy  dz 

z=0y=0x=0 

Then  from  Eqns.  (40)  and  (41),  we  can  compute  the  macroscale  CTE’s  as: 

(/I 


(41) 


(42) 


where  {/}  is  given  by  the  expression: 


c b a 

{/}  = I I ^[C\[a}  dx  dy  dz 

z=0y=Qx=0 


(43) 


3.2  Continuum  Results  using  SAM 

A code  called /^TE';i'-20  (pronounced  as  microtech)  was  written  in  FORTRAN  77  to 
implement  SAM.  The  code  was  executed  to  estimate  the  thermo-elastic  constants  for  the 
seven  examples,  whose  constituent  material  properties  are  listed  in  Table  1 . Input  to  the  code 
were  the  unit-cell  dimensions,  yarn  geometry  information,  constituent  material  properties. 
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and  the  number  of  divisions  required  to  discretize  the  unit-cell  in  the  x,  and  z directions. 
The  element  stiffness  matrix  [C]  was  determined  by  computing  the  elasticity  matrix  for  the 
material  point  at  the  geometric  center  of  the  element,  and  transforming  it  to  the  unit-cell 
coordinate  system.  The  predicted  macroscale  stiffness  matrix  [C^],  and  consequently  the 
macroscale  compliance  matrix  will  not  be  symmetric  due  to  the  approximate  nature  of  the 
analysis.  Therefore  the  macroscale  compliance  matrix  was  made  symmetric  by  averaging 
the  off-diagonal  compliance  coefficients.  The  macroscale  elastic  constants  were  computed 
by  comparing  the  symmetrized  compliance  coefficients  with  that  of  a homogenous, 
orthotropic  medium. 

The  results  for  a bimaterial  medium  (example  2)  are  given  in  Table  15.  The  results 
for  example  1 are  not  listed,  since  it  is  obvious  that  SAM  would  predict  the  elastic  constants 
for  an  isotropic  medium  exactly.  The  bimaterial  medium  consisted  of  two  different  layers 
of  isotropic  materials  of  equal  thickness  altematingly  stacked  in  the  z-direction.  The  elastic 
constants  for  the  bimaterial  medium  were  compared  with  an  exact  solution  (derivation 
explained  in  Section  2.1.3),  and  with  the  previously  computed  finite  element  results.  It  can 
be  observed  that  SAM  marginally  under-predicts  the  longitudinal  and  transverse  Young’s 
moduli,  while  the  inplane  and  transverse  shear  moduli  are  exact. 


Table  15.  Continuum  properties  for  example  2 using  SAM. 


Ex,Ey 

(GPa) 

Ez 

(GPa) 

Gxz, 

Gyz 

(GPa) 

Cry 

(GPa) 

^yz 

V;cy 

X 

10-6/oC 

Oz^ 

X 

10-6/oC 

Example  2 
(bimaterial 
medium) 

fxTEx-20 

(SAM) 

36.02 

8.72 

2.48 

15.23 

0.599 

0.183 

3.88 

52.20 

^iTEx-lO 

(FEA) 

36.79 

9.79 

2.48 

15.23 

0.312 

0.208 

8.19 

59.60 

exact 

solution 

36.79 

9.79 

2.48 

15.23 

0.312 

0.208 

8.19 

59.60 
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Table  16  presents  the  SAM  results  for  two  cases  of  unidirectional  composite 
(examples  3 and  4).  The  fiber  and  matrix  had  identical  Poisson’s  ratio  in  example  3,  and 
different  Poisson’s  ratio  in  example  4.  The  SAM  results  were  compared  with  the  finite 
element  results  from  the  previous  chapter,  and  with  analytical  solutions  for  unidirectional 
composite  properties.  The  analytical  expressions  used  were  the  rule  of  mixtures  formulae 
for  El  and  vn  and  the  Halpin-Tsai  equations  (Halpin  and  Tsai,  1969)  for  Ej,  Gn  and  vjt. 
All  of  the  unidirectional  composite  thermo-elastic  constants  but  for  Ej  and  aj  were  found 
to  match  well  with  the  compared  data.  Table  17  compares  the  SAM  results  for  three  textile 
composites  (examples  5, 6 and  7)  with  available  results.  In  all  three  cases  the  thermo-elastic 
constants  obtained  by  implementing  SAM  were  in  good  agreement  with  the  available  results. 


Table  16.  Continuum  properties  for  examples  3 and  4 using  SAM. 


El 

(GPa) 

Ej 

(GPa) 

Glt 

(GPa) 

Gjt 

(GPa) 

vlt 

vtt 

at 

X 

10-6  PC 

OLj 

xlO-6/°C 

fiTEx-20 

(SAM) 

64 

50.24 

10.45 

8.36 

0.341 

0.300 

12.41 

21.51 

Example  3 

(unidirect. 

composite) 

fiTEx-lO 

(FEA) 

63.55 

36.48 

12.93 

9.94 

0.300 

0.232 

15.74 

40.79 

rule  of 
mixt./ 
Halpin- 
Tsai  eqns. 

64 

34.55 

11.26 

13.29 

0.300 

0.300 

15.63 

55.11 

fj.TEx-20 

(SAM) 

43.35 

32.47 

4.13 

3.04 

0.245 

0.218 

7.40 

11.60 

Example  4 

(unidirect. 

composite) 

fiTEx-\0 

(FEA) 

43.12 

18.15 

5.59 

3.92 

0.242 

0.222 

7.40 

25.44 

rule  of 
mixt./ 
Halpin- 
Tsai  eqns. 

43.40 

14.79 

4.45 

5.91 

0.260 

0.252 

6.77 

34.24 
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Table  17.  Continuum  properties  for  examples  5,  6 and  7 using  SAM. 


(GPa) 

Ez 

(GPa) 

Gxz, 

Gyz 

(GPa) 

(GPa) 

Vxz^ 

Vyz 

Vxy 

ax‘^,tty‘^\ 

10-6/oc 

Oz^\ 

10-6/oC 

Example  5 

(plain- 

weave) 

fiTEx-2Q 

(SAM) 

12.46 

6.62 

1.64 

1.67 

0.399 

0.162 

29.10 

68.48 

fiTEx-lO 

(FEA) 

11.81 

6.14 

1.84 

2.15 

0.408 

0.181 

28.36 

79.57 

Dasgupta 

results 

14.38 

6.25 

1.94 

3.94 

0.463 

0.167 

22.50 

86.00 

Example  6 

(plain- 

weave) 

fiTEx-20 

(SAM) 

63.41 

11.13 

3.79 

4.24 

0.402 

0.027 

1.36 

21.53 

fiTEx~\Q 

(FEA) 

53.61 

10.88 

4.41 

4.72 

0.365 

0.128 

1.56 

22.71 

TEXCAD 

64.38 

11.49 

5.64 

4.87 

0.396 

0.027 

1.33 

20.71 

Example  7 

(5-harness 

weave) 

fiTEx-20 

(SAM) 

69.30 

11.62 

4.06 

4.73 

0.355 

0.031 

1.21 

20.25 

fiTEx-lO 

(FEA) 

64.51 

11.33 

4.45 

4.85 

0.329 

0.047 

1.55 

22.03 

TEXCAD 

66.33 

11.51 

4.93 

4.89 

0.342 

0.034 

1.46 

21.24 

3.3  Selective  Averaging  Method  (SAMI  for  Plate  Properties 

In  this  section,  the  SAM  procedure  to  compute  the  plate  stiffness  coefficients  and 
plate  thermal  coefficients  (for  a thin  textile  composite)  is  described.  To  distinguish  between 
the  macrolevel,  mesolevel  and  microlevel  [A],  [5]  and  [D]  matrices,  an  over-tilde  is  used  to 
denote  the  mesolevel  stiffness,  and  a superscript  ”M”  is  used  to  denote  the  macrolevel  plate 
stiffness.  However  in  the  remaining  sections,  [A],  [fi]  and  [D]  (without  a superscript  or 
overscript)  will  refer  to  the  macroscale  plate  stiffness  matrices.  Also  in  this  section,  the 
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complete  plate  stiffness  matrix  on  the  macroscale,  as  defined  by  Eqn.  (23)  will  be  denoted 
by  [C^],  such  that: 


[C^] 


(44) 


The  procedure  to  find  the  plate  stiffness  and  thermal  coefficients  is  analogous  to  that 
used  to  find  the  continuum  thermo-elastic  constants.  To  find  the  first  column  of  the  effective 
plate  stiffness  matrix,  the  unit-cell  is  discretized  into  slices  (mesolevel)  and  elements 
(microlevel)  as  shown  in  Fig.  20.  The  unit-cell  is  subject  to  the  deformation  given  by  e^=l- 
The  following  assumptions  are  made  regarding  the  midplane  strains  and  curvatures  on  the 
macrolevel,  mesolevel  and  microlevel: 


£«  = = 0 i = 2,3 

xf  = = x^  = Q / = 1,3 


(45) 


It  is  also  assumed  that  the  non-zero  strain  component  e^o,  and  the  force  resultant  Nx  are 
uniform  within  the  mesoscale  and  macroscale  respectively.  These  two  assumptions  can  be 
expressed  as  the  following  equations: 


^;c0  “ 

Nx  = 


(46) 


The  mesolevel  stiffness  coefficient  A j j can  then  be  obtained  by  averaging  the  corresponding 
element  stiffness  coefficients  over  the  slice  (consequence  of  the  isostrain  assumption)  as: 

c b 

^ i J j 2ii(Jc,y,z)  dy  dz  (47) 

z—0y=0 

where  i is  the  plane-stress  stiffness  coefficient  in  the  classical  lamination  theory  (Agarwal 
and  Broutman,  1990),  which  has  been  transformed  to  the  unit-cell’s  xyz-coordinates  (for  an 
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isotropic  material,  = — — ).  The  macroscale  force  and  moment  resultants  can  be 

(1  - v2) 

expressed  in  terms  of  the  microscale  stresses  by  the  following  relations: 


c b a 


= -!- 


= -V 


^ / 1 f"' 

z = 0 y = 0j:  = 0 
c b a 

I 

z=0y=0x=0 


dx  dy  dz 


i = 1,2,3 


(48) 


dx  dy  dz 


The  assumption  of  uniform  force  resultants  on  the  macroscale  and  Eqn.  (48)  yields  the 
following  expressions  for  the  first  column  of  plate  stiffness  coefficients: 


a 


(49) 


1 

ab 


c b a 


III 

z=0y=0x=0 


jML 

iii(x) 


Qn(x,y,z)  dx  dy  dz 


a = 1,2,3) 


(50) 


(~'M 

Si 


Qi\{x,y,z)  dx  dy  dz 


z=0y=0x=0 


(i  = 1,2,3) 
7 = 1 + 3 


(51) 


A similar  procedure  is  followed  to  compute  the  second  column  of  the  effective  plate  stiffness 
matrix. 


In  the  case  of  shear  loading  = 1,  the  unit-cell  is  discretized  into  slices  parallel 
to  the  yz-plane.  The  assumptions  of  isostrain  and  uniform  force  resultants  are  reversed  from 
the  case  of  normal  loading.  The  force  resultant  Nxy  is  assumed  to  be  uniform  within  a slice, 
such  that 


= Nxy 


(52) 
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It  is  also  assumed  that  jxy  is  the  only  non-zero  deformation  component  on  the  macrolevel, 
mesolevel  and  microlevel.  Averaging  the  element  compliance  coefficients  over  the  slice,  we 
obtain: 


1 


b a 

ib  j j 

y=0x=0 


1 

Q23(x,y,z) 


dx  dy 


(53) 


The  mesolevel  stiffnesses  are  then  averaged  over  the  volume  of  the  unit-cell  to  yield  the  third 
column  of  the  plate  stiffness  matrix,  as  follows 

c b a 


cM  _ 1 

ab 


, , , 


{i  = 1,2,3) 


(54) 


z = 0>"  = 0jc=0 


c b a 

rM  = ± [ [ [ , Qi3(^,y,z) 

ab  j j j Q33(x,y,z) 

z=0y=0x=0 


Q33(z)  dx  dy  dz 


(i  = 1,2,3) 
7 = 1 + 3 


(55) 


The  expressions  for  the  fourth,  fifth  and  sixth  columns  of  the  plate  stiffness  matrix 
can  be  obtained  in  a similar  fashion  to  that  explained  above,  except  that  the  one  of  the 
curvature  component  will  be  non-zero  instead  of  a midplane  strain  component.  For  example, 
to  determine  the  fourth  column,  Xx  will  be  the  only  non-zero  deformation  on  the  macroscale, 
mesoscale  and  microscale.  Assuming  that  the  curvature  is  uniform  within  a slice,  we  get 
Hx  = Hx-  Then  by  averaging  the  element  stiffness  coefficients  over  the  slice,  we  obtain  an 
expression  analogous  to  Eqn.  (47)  as: 


Dii(x) 


c b 

ill 

z = 0>’  = 0 


Z^Qu(x,y,z)  dy  dz 


(56) 


The  moment  resultant  Mx  is  assumed  to  be  uniform  on  the  mesoscale  such  that  = Mx- 
By  averaging  the  slice  compliance  coefficients  over  the  unit-cell  volume,  we  get  the 
following  relations  for  the  fourth  column  of  the  plate  stiffness  matrix: 
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a 


(57) 


^i4 


J_ 

ab 


c b a 

IK 

z = 0y  = 0A:  = 0 


dm 

11 

D^^{x) 


Qi\(x,y,z)  dx  dy  dz  (i  = 1,2,3) 


(58) 


/~<M  1 

K ab 


c b a 


HI- 

z=0y=0x=0 


Dii(x) 


Qi^(x,y,z)  dx  dy  dz 


a = 1,2,3) 
7 = 1 + 3 


(59) 


The  fifth  and  sixth  columns  of  the  plate  stiffness  matrix  can  be  found  by  using  a similar 
procedure. 

To  find  the  plate  thermal  coefficients,  and  are  assumed  as  zero  in  the 
macrolevel,  mesolevel  and  microlevel.  The  thermal  stresses  developed  in  the  microscale  due 
to  a uniform  temperature  difference  of  A T=Tq  are  given  by: 


{^7}  = - [Q]{a}Tf, 


(60) 


where 


Ox' 

'ax' 

{a}  = ■ 

Oy 

^xy 
» .. 

{a}  = ■ 

ay 

axy 

The  macroscale  plate  constitutive  equation  will  reduce  to 


(61) 


(62) 


By  averaging  the  microscale  stresses  given  by  Eqn.  (60)  over  the  unit-cell  volume  (using 
Eqn.  48),  and  equating  to  the  macroscale  force  and  moment  resultants  in  Eqn.  (62),  we  get 
the  following  relations  for  the  plate  CTE’s: 


h 


(63) 
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where  I\  and  I2  are  integrals  given  by  the  expressions: 


c b a 


III 


[Q]{a}  dx  dy  dz 


z = 0y  = 0;c  = 0 


c b a 


(64) 


III 


z [Q]{a}  dx  dy  dz 


z = 0>»  = 0x=0 


3.4  Plate  Results  using  SAM 


The  plate  coeffieients  and  plate  CTE’s  were  computed  for  examples  2-7  by 
implementing  SAM,  and  compared  with  the  finite  element  results  presented  in  the  previous 
chapter.  The  computed  coefficients  for  the  examples  are  listed  in  Tables.  18-20.  The 
properties  for  the  bimaterial  plate  were  predicted  very  accurately  (Table  18).  For  the 
unidirectional  composite  examples,  SAM  was  able  to  predict  the  [A]  matrix  coefficients 
(except  for  A66),  D\\^  and  the  plate  CTE’s  very  well  (Table  19).  However  SAM  grossly 
over-predicted  D\2  and  D22  and  under-predicted  the  coefficient  D66-  For  the  textile 
composite  examples  (Table  20)  SAM  predicted  all  but  A 12  andDi2  with  very  good  accuracy. 
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Table  18.  Non-zero  [A],  [5]  and  [D]  coefficients  for  example  2 (bimaterial  plate)  using 
SAM. 


^1I>  ^22 
xl06 

A/2 

x1Q6 

A56 

xl06 

Bii,  B22 
xl03 

B]2 

xl03 

B66 

xl0“3 

fiTEx-20 

(SAM) 

9.844 

2.045 

3.899 

- 0.565 

-0.108 

- 0.228 

fiTEx-lO 

(FEA) 

9.832 

2.043 

3.895 

- 0.563 

-0.108 

- 0.228 

lamination  theory 
for  two  plies 

9.832 

2.043 

3.895 

- 0.563 

-0.108 

- 0.228 

lamination  theory 
using  3-D  elastic 
constants 

9.844 

2.048 

3.899 

0 

0 

0 

B>11,  E>22 
xlO-3 

D]2 

xlO-3 

B>66 

xl0“^ 

a^P,  ayP 

xlO-6/oc 

PxP.^yP 

/«C/m 

fiTEx-20 

(SAM) 

53.727 

11.163 

21.282 

17.828 

0.170 

f^TEx-lO 

(FEA) 

53.590 

11.149 

21.220 

17.800 

0.170 

lamination  theory 
for  two  plies 

53.573 

11.131 

21.220 

17.814 

0.170 

lamination  theory 
using  3-D  elastic 
constants 

53.762 

11.183 

21.293 

8.190 

0 

Note:  [A],  [5]  and  [D]  coefficients  in  SI  units 
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Table  19.  Non-zero  [A],  [5]  and  [D]  coefficients  for  examples  3 and  4 (unidirectional 
composite)  using  SAM. 


An 

x1Q6 

An 

xl06 

A22 

xl06 

A66 

xl06 

Example  3 

uTEy-20 

(SAM) 

0.688 

0.175 

0.475 

0.102 

fiTEx-lO 

(FEA) 

0.690 

0.149 

0.496 

0.177 

Halpin-Tsai  Eqns. 
and  lamination 
theory 

0.673 

0.109 

0.363 

0.113 

Example  4 

fiTEx-20 

(SAM) 

0.443 

0.074 

0.261 

0.040 

fiTEx-lO 

(FEA) 

0.452 

0.062 

0.285 

0.114 

Halpin-Tsai  Eqns. 
and  lamination 
theory 

0.444 

0.039 

0.151 

0.045 

Dn 

xlO-6 

Di2 

xlO-6 

D22 

xlO-6 

D66 

xlO-6 

ttxPx 

10-6 

ttyPX 

10-6 

Example  3 

fiTEx~20 

(SAM) 

3.750 

1.029 

3.112 

0.541 

14.476 

24.750 

luTEx-\0 

(FEA) 

3.589 

0.596 

1.980 

0.947 

15.489 

26.184 

Halpin-Tsai 
Eqns.  and 
lamination 
theory 

5.606 

0.908 

3.026 

0.939 

15.625 

55.112 

Example  4 

uTEx-20 

(SAM) 

2.308 

0.446 

1.799 

0.195 

6.628 

13.076 

f^TEx-W 

(FEA) 

2.256 

0.224 

0.873 

0.568 

7.378 

13.188 

Halpin-Tsai 
Eqns.  and 
lamination 
theory 

3.702 

0.328 

1.262 

0.371 

6.774 

34.239 

Note;  [A],  [fi]  and  [D]  coefficients  in  SI  units 
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Table  20.  Non-zero  [A],  [6]  and  [D]  coefficients  for  examples  5,  6 and  7 using  SAM. 


An , A22 

X 106 

An 

X 106 

A66 

X 106 

Bn 

xl03 

Example  5 

HTEx-20  (SAM) 

2.667 

0.446 

0.379 

0 

fiTEx-\0  (FEA) 

2.681 

0.565 

0.489 

0 

lamination  theory 
using  3-D  constants 

2.783 

0.503 

0.490 

0 

Example  6 

HTEx-2Q  (SAM) 

12.215 

0.577 

1.095 

0 

HTEx-\0{¥EA) 

12.090 

3.470 

1.223 

0 

lamination  theory 
using  3-D  constants 

13.938 

1.787 

1.208 

0 

Example  7 

^TEx-20  (SAM) 

16.039 

0.631 

1.209 

0.590* 

fiTEx-lO  (FEA) 

14.683 

1.351 

1.210 

0.495* 

lamination  theory 
using  3-D  constants 

16.531 

0.770 

1.239 

0 

Dji,  E>22 

X 10-3 

Di2 

X 10-3 

Dee 

X 10-3 

GxP,  GyP 
X 10-6/oC 

PxP 

Pdm 

Example  5 

liTEx-20  (SAM) 

emi 

1.590 

1.360 

27.505 

0 

fxTEx-\0  (FEA) 

5.687 

1.518 

1.577 

27.465 

0 

lamination  theory 
using  3-D 
constants 

12.054 

2 All 

2.124 

28.363 

0 

Example  6 

juTEx-20  (SAM) 

44.782 

2.733 

4.398 

1.984 

0 

juTEx-lO  (FEA) 

41.695 

0.373 

5.879 

1.480 

0 

lamination  theory 
using  3-D 
constants 

75.942 

9.734 

6.582 

1.556 

0 

Example  7 

fiTEx-20  (SAM) 

81.164 

3.301 

5.868 

2.422 

-0.028* 

^TEx-\0(FEA) 

90.072 

1.123 

6.149 

2.910 

-0.037* 

lamination  theory 
using  3-D 
constants 

87.283 

4.195 

6.753 

1.550 

0 

* In  example  7,  £22=  -^n  and 

Note:  [A],  [B]  and  \D]  coefficients  in  SI  units 


CHAPTER  4 

FINITE  ELEMENT  MODELS  FOR  STRENGTH  PROPERTIES 

In  Chapter  2,  we  had  demonstrated  finite  element  procedures  to  model  a general 
textile  composite  either  as  a three-dimensional  (continuum)  material  or  as  a thin  plate/beam 
to  predict  their  corresponding  thermo-mechanical  coefficients.  In  this  chapter,  we  extend  the 
same  numerical  models  to  compute  the  thermal  residual  stresses  due  to  processing  in  the 
yarns  and  the  matrix.  Then  the  numerical  models  are  used  to  study  the  strength  behavior  of 
the  composite  by  predicting  the  failure  envelopes  for  thin  and  thick  textile  composites. 

4. 1 Thermally  Induced  Residual  Microstresses 

The  thermal  residual  microstresses  are  induced  in  the  yarn  and  matrix  materials  due 
to  the  mismatch  in  their  corresponding  CTE’s.  The  difference  between  the  composite  curing 
temperature  and  room  temperature  then  serves  as  the  driving  force  to  create  these 
microstresses.  Since  composites  designed  for  high  temperature  applications  are  fabricated 
at  higher  temperatures,  the  residual  microstresses  become  relevant  in  the  strength 
considerations  of  such  composites.  The  residual  microstresses  in  the  vicinity  of  the 
yarn-matrix  interface  are  particularly  important  as  they  could  lead  to  failure  due  to 
debonding. 

Determination  of  residual  microstresses.  Let  To  be  the  difference  between  room 
temperature  and  the  composite  fabrication  temperature.  Since  the  composite  is  stress  free  at 
the  fabrication  temperature,  which  is  above  room  temperature,  Tg  is  generally  negative.  The 
residual  microstresses  in  the  yarn  and  the  matrix  are  obtained  by  superposing  the 
microstresses  due  to  the  two  load  cases  as  explained  below.  In  the  first  load  case,  the  unit-cell 
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is  constrained  from  expanding  by  fixing  the  comer  nodes  of  the  unit-cell  and  enforcing  zero 
displacement  difference  between  corresponding  nodes  on  opposite  faces  of  the  unit-cell 
(periodic  displacement  EC’s).  A temperature  difference  Tg  is  applied  to  all  elements  in  the 
finite  element  model.  This  is  exactly  the  same  problem  we  solve  for  finding  the 
three-dimensional  (continuum)  CTE’s.  The  applied  boundary  conditions  mean  that  all  the 
macrostrain  components  are  equal  to  zero  = 0 ,AT’^  = Tq).  Then  the  corresponding 

macrostresses  required  to  restrain  the  unit-cell  expansion  are  given  by: 

K)  = - [C]{a^)To  (65) 

In  the  second  load  case,  deformations  are  applied  so  as  to  reverse  the  macrostresses 
developed  in  the  first  load  case.  This  can  be  accomplished  by  imposing  the  deformations 
{e^}  = {a^]  Tg  and/d  r^=0.  It  can  be  noted  that  the  macrostresses  developed  in  the  second 
loading  case,  given  by  [C]{a^}7^  are  equal  and  opposite  to  the  macrostresses  in  Eqn.  (65). 
The  microstresses  from  both  load  cases  are  superposed  to  obtain  the  residual  stresses  due  to 
free  thermal  expansion. 

The  same  idea  can  be  extended  to  finding  the  residual  microstresses  in  the  plate 
model.  Then  the  deformations  to  be  applied  in  the  first  load  case  are  { e^)  =0,{x^}=0  and 
AT^  =Tg  \ and  the  deformations  in  the  second  load  case  are  {Eq}  = {aP}  Tg  , {x^}  = {^p} 
Tg  andzl  = 0.  The  residual  microstresses  were  computed  for  the  plain-weave  textile  beam 
at  the  Gaussian  center  of  the  elements  in  the  unit-cell.  The  beam  was  assumed  to  be  in  the 
x-z  plane  with  unit-cells  repeating  in  the  x-direction.  figure  21  shows  the  thermal  stress 
contours  for  Oxx,  o^zz  Txz-  The  composite  curing  temperature  was  assumed  to  be  150°C 
above  room  temperature. 
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(a) 


(b) 


(c) 


Figure  21.  Thermal  stress  contours  in  a plain-weave  beam  for  zl  7'=-150°C. 
(a)Oxr;  (b)o-2,;  (c)r^,. 
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4.2  Strength  Modeling  of  Textile  Composites 

There  are  many  failure  criteria  or  strength  theories  for  unidirectional  fiber 
composites.  This,  for  example,  includes  maximum  stress  theory,  maximum  strain  theory  and 
Tsai-Hill  theory  (Agarwal  and  Broutman,  1990).  Even  though  failure  of  a material  is  a very 
complex  phenomenon,  engineering  strength  theories  such  as  mentioned  above  have  been 
found  to  be  useful  in  design.  The  interpretation  of  strength  values  obtained  from  such 
theories  are  different  for  different  materials.  For  example  in  metal  matrix  composites  the 
failure  envelope  obtained  using  the  above  theories  will  correspond  to  the  initial  yield  surface 
(Dvorak  et  al.,  1973).  In  graphite/epoxy  composites  the  failure  theories  can  be  used  to  predict 
fiber  or  matrix  failure.  In  the  present  study  our  intent  is  to  explore  the  possibility  of 
developing  such  failure  criteria  for  textile  composites. 

4.2.1  Determination  of  Composite  Failure  Envelope 

Our  approach  is  similar  to  that  used  by  Dvorak  et  al.  (1973).  A state  of  homogeneous 
deformation,  corresponding  to  each  of  the  six  macrostrain  components,  are  independently 
applied  to  the  unit-cell  by  imposing  the  boundary  conditions  explained  in  Section  2. 1 . For 
each  case,  the  various  stress  components  are  computed  in  the  elements  in  the 
unit-cell — typically  at  the  element  Gaussian  integration  points.  These  stresses  will  be 
referred  to  as  microstresses.  Assuming  linear  elastic  behavior,  the  microstresses  can  be 
computed  for  any  arbitrary  combination  of  the  macrostrain  components.  Since  we  know  the 
macroscale  elasticity  matrix,  we  can  find  a relation  between  microstresses  at  a point  and  any 
arbitrary  state  of  macrostress  as: 

{o-}  = [F]  [o^]  (66) 

[F]  can  be  considered  as  a matrix  of  influence  coefficients,  which  is  evaluated  at  the 
integration  points  of  all  the  elements  in  the  unit-cell.  We  also  assume  that  the  failure  behavior 
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of  the  matrix  material  and  the  yam  is  known.  For  instance,  let  the  failure  criterion  of  the 
matrix  be  of  the  form  [H]  {<7}  matrix  = 1-  Then  the  failure  criterion  for  the  composite  is 
obtained  from  Eqn.  (66)  as  [H]  [F]  {a^}  = 1.  The  same  idea  also  applies  for  the  yarn.  The 
textile  composite  is  assumed  to  have  failed  if  there  is  failure  on  the  microscale  in  any  one 
of  the  constituent  materials — either  matrix  or  the  yam.  By  varying  the  macrostresses  using 
a numerical  simulation,  failure  envelopes  can  be  obtained  for  the  idealized  homogeneous 
material  (Fig.  22a).  It  might  be  noted  that  Eqn.  (66)  can  be  modified  to  include  the  thermal 
residual  stress  field  in  the  unit-cell  as 

[a]  = [F]  [a^]  - \ar\To  (67) 

where  {07}  is  the  matrix  of  thermal  microstresses  computed  at  the  element  integration  point 
for  a unit  temperature  difference. 


Figure  22.  Composite  failure  envelopes. 

(a)  continuum  failure  envelope  in  the  space  of  macrostresses; 

(b)  plate  failure  envelope  in  the  space  of  force  and  moment  resultants. 
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4.2.2  Effect  of  Stress  Gradients  on  Strength 


The  strength  analysis  for  a three-dimensional  composite  can  be  extended  for  thin 
composites  using  the  plate  model.  As  mentioned  in  Section  2.2,  the  macrostresses  will  not 
be  homogeneous  thorough  the  thickness  in  such  composites.  Then  the  composite  failure  will 
be  determined  by  the  stress  gradients  through  the  thickness,  represented  by  the  averaged 
force  resultants  (AO  and  moment  resultants  (A/).  The  composite  failure  criterion  for  the  plate 
model  will  be  of  the  form: 


mm 


N 

M 


= 1 


(68) 


Thus  the  failure  envelope  of  the  composite  will  be  in  the  six-dimensional  space  of  the  force 
resultants  and  the  moment  resultants  (Fig.  22b).  The  above  procedure  was  demonstrated 
using  the  beam  model,  for  the  case  of  a plain-weave  textile  composite. 


Failure  envelope  for  textile  composite  beam.  For  a textile  composite  beam,  the 
failure  envelope  is  constructed  in  the  space  of  the  three  force  resultants  P,  M and  V.  The  textile 
composite  beam  is  assumed  to  be  in  the  jcz-plane  with  unit-cells  repeating  in  the  x-direction. 
The  unit-cell  is  discretized  with  eight-node  isoparametric  plane  strain  finite  elements.  Three 
linearly  independent  deformations,  as  explained  in  Section  2.3,  are  applied  to  the  unit-cell: 
(a)  unit  axial  strain;  (b)  unit  curvature  such  that  transverse  shear  strain  is  zero;  (c)  unit 
transverse  shear  strain.  For  each  deformation,  from  the  finite  element  results,  the 
microstresses  cr^,  and  txz  are  computed  at  the  Gaussian  center  of  each  element.  The 
microstresses  in  the  element  for  a combination  of  loads  (deformations)  are  given  by: 


r„M'i 

< (^‘zz 

> 

= [A']  < 

\ 

y? 

(69) 
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where  the  each  column  of  [A']  corresponds  to  the  microstresses  for  unit  deformation.  For 
example,  the  first  column  of  [A']  represents  the  three  microstress  components  in  the 
element  for  unit  Eq.  We  also  know  that 


where  [S]  is  the  beam  compliance  matrix.  From  Eqns.  (69)  and  (70)  we  arrive  at  an 
expression  analogous  to  Eqn.  (66): 


where  [F]  = [A']  [5].  The  composite  is  assumed  to  have  failed  if  there  is  failure  in  any  one 
of  the  finite  elements  in  the  unit-cell.  Matrix  failure  is  determined  by  either  maximum 
principal  stress  or  von  Mises  criteria.  Since  the  yarn  is  assumed  to  be  transversely  isotropic, 
yarn  failure  is  determined  by  either  maximum  strain  theory  for  a unidirectional  composite 
or  the  Tsai-Wu  criteria. 

Beam  Failure  Envelope  Results.  The  strength  properties  used  for  the  constituent 
materials  in  the  beam  are  as  follows; 

Yarn:  ^=1725  MPa,  Oi  1 366  MPa,  Ot'^=A2  MPa,  aT^=23Q  MPa,  Xit=95  MPa 

Matrix:  cr^=70  MPa,  a^=100  MPa 

where  the  superscripts  ’ T and  ’ C’  refer  to  the  tensile  and  compressive  strengths  respectively. 
The  failure  envelopes  were  developed  using  two  different  sets  of  failure  criteria.  In  the  first 
case  maximum  principal  stress  criterion  was  used  to  determine  the  matrix  failure  and  the 
maximum  strain  theory  for  unidirectional  fiber  composite  was  used  for  the  yam.  In  the 
second  case  von  Mises  criterion  was  used  for  the  matrix  and  the  Tsai-Wu  criteria  for  the  yam. 
Both  structural  (beam)  and  continuum  failure  envelopes  were  developed.  As  was  explained 


P' 

[S]  \m  . 
V 


(70) 
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earlier,  in  the  beam  model,  periodic  EC’s  were  imposed  between  corresponding  nodes  on 
the  left  and  right  ends  of  the  unit-cell.  In  the  continuum  model,  the  periodic  EC’s  were  in 
addition  imposed  between  the  top  and  bottom  surfaces.  The  failure  envelopes  in  Figures  (23) 
through  (26)  were  obtained  using  the  first  set  of  failure  criteria.  Figures  (23)  and  (24)  depict 
the  structural  failure  envelope  in  the  P-M  space  based  on  yarn  and  matrix  failure 
respectively.  Figures  (25)  and  (26)  illustrate  similar  continuum  failure  envelopes  in  the  space 
of  the  macroscale  normal  stresses  in  the  x-  and  z-directions.  As  expected,  the  envelope 
reduced  in  size  with  increasing  shear  force  resultant.  If  we  assume  that  the  beam  is  made  of 
a homogeneous  but  orthotropic  material  with  properties  as  predicted  by  the  continuum 
model,  then  one  can  derive  structural  failure  envelopes  from  the  continuum  failure  envelopes 
using  simple  beam  theories.  The  derived  structural  failure  envelopes  are  compared  with  that 
obtained  from  direct  micromechanical  analyses  in  Figures  (27)  and  (28).  One  can  note  that 
the  continuum  failure  criteria  are  very  conservative  for  the  case  of  a thin  beam.  Figures  (29) 
and  (30)  show  similar  comparisons  for  the  second  set  of  results  obtained  using  the  quadratic 
failure  criteria  for  the  constituent  materials. 
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Figure  23.  Beam  failure  envelope  based  on  yam  failure,  (maximum  strain  theory) 


150 


100  - 


50  - 


c 0 - 


• . o «• 


V=0 

OOOOO  V =5x10* 

v=lxl0^ 


-50  L 


-100- 


-150 


Figure  24.  Beam  failure  envelope  based  on  matrix  failure,  (maximum  principal 
stress  criteria) 


73 


N 

S 

to 


xlO> 

0.5 1 

0 - 

-0.5  - 

-1  - 

-1.5  - 

-2  - 

-2.5  - 

-V 

-4  -3  -2  -1  0 


= 0 

00000  r„=25  MPa 
t=50  MPa 


. O* 
.0* 
.0* 
.0* 

.0* 

.0* 


1 


a 


\x 


Xl08 


Figure  25.  Continuum  failure  envelope  based  on  yarn  failure,  (maximum  strain 
theory) 
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Figure  26.  Continuum  failure  envelope  based  on  matrix  failure.  ( maximum  principal 
stress  criteria) 
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Figure  27.  Comparison  of  failure  envelopes  based  on  yam  failure,  (maximum  strain 
theory) 


a.xial  force  (P) 

Figure  28.  Comparison  of  failure  envelopes  based  on  matrix  failure,  (maximum 
principal  stress  criteria) 
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Figure  29.  Comparison  of  failure  envelopes  based  on  yarn  failure.  (Tsai-Wu  criteria) 


Figure  30.  Comparison  of  failure  envelopes  based  on  matrix  failure.  (von-Mises 
criteria) 


CHAPTER  5 

ISSUES  IN  MICROMECHANICAL  MODELING 


The  complex  yarn  architectures  in  textile  composites  make  traditional  finite  element 
modeling  very  difficult.  Traditional  finite  element  models  assume  that  the  material 
properties  are  constant  or  vary  smoothly  within  an  element.  In  the  context  of  textile 
composites,  it  means  that  the  yam  and  matrix  materials  are  modeled  by  separate  elements, 
with  common  nodes  at  the  yarn-matrix  interface.  This  is  indeed  preferable  because  the 
stresses  at  the  interface  can  be  computed  accurately.  However  meshing  the  interstitial  matrix 
region  becomes  time-consuming  and  difficult.  This  region  which  is  essentially  a collection 
of  multiply-connected  matrix  pockets,  requires  a very  fine  mesh  to  capture  the  pocket 
geometry.  Dasgupta  et  al.  (1990)  and  Whitcomb  (1991)  presented  detailed  finite  element 
analyses  of  the  unit-cell  for  plain-weave  architectures.  In  this  chapter,  alternative  finite 
element  modeling  techniques  (Marrey  and  Sankar,  1994)  are  discussed,  which  are  valid  for 
any  yam  architecture.  The  methods  were  tested  with  simple  two-dimensional  examples,  and 
compared  with  the  displacement  fields  using  traditional  finite  elements.  Some  aspects  of 
finite  element  mesh  generation,  in  the  context  of  textile  composites,  are  also  discussed. 

5.1  Finite  Element  Modeling  of  the  Unit-Cell 

Traditional  finite  elements  use  homogeneous  elements,  i.e.,  elements  that  are 
comprised  of  only  one  material.  Two  methods  to  circumvent  traditional  finite  element 
modeling  are  described  in  the  following  section.  The  first  uses  incompatible  elements,  which 
employs  homogeneous  elements  for  the  yarn  and  matrix,  though  there  is  a node  mismatch 
at  the  yarn-matrix  interface.  The  second  method  utilizes  inhomogeneous  elements  to  model 
the  unit-cell. 
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5.1.1  Incompatible  Elements 

The  following  method  is  a modification  of  previously  published  work  on 
independently  modeling  substructures  with  finite  elements  (Ransom  et  al.,  1993).  The 
matrix  and  each  of  the  yarns  are  meshed  independently  with  finite  elements.  In  general,  the 
nodes  on  the  surface  of  the  yarns  will  not  coincide  with  the  matrix  nodes  on  the  yarn-matrix 
interface.  In  fact,  most  textile  composite  unit-cells  due  to  multidirectional  reinforcements, 
contain  multiply-connected  matrix  pockets.  These  matrix  regions  require  an  extremely  fine 
mesh  to  capture  their  geometry,  whereas  the  yarns  can  be  modeled  with  a coarse  mesh.  This 
method  has  therefore  two  advantages — the  first  being  the  relative  ease  of  generating 
independent  meshes  for  the  yams  and  the  matrix.  Secondly,  the  effective  degrees  of  freedom 
in  the  numerical  model  are  reduced,  since  the  yarns  are  modeled  with  a coarse  mesh. 

Figure  31(a)  shows  an  example  of  a rectangular  unit-cell  with  an  ellipsoidal 
inclusion.  The  matrix  mesh  is  shown  in  Fig.  31(b)  and  the  inclusion  mesh  is  shown  in  Fig. 
31(c).  The  matrix  degrees  of  freedom  on  the  inclusion-matrix  interface  are  denoted  as  qh'”; 
the  matrix  degrees  of  freedom  in  the  interior  as  and  the  inclusion  degrees  of  freedom  are 
denoted  as  qy.  The  total  strain  energy  in  the  unit-cell  is  given  by 


where  Ky  is  the  finite  element  stiffness  for  the  inclusion;  and  is  the  stiffness  for  the  matrix 
material,  which  is  divided  into  four  sub-matrices.  Ry  and  are  the  external  load  vectors  in 

the  inclusion  and  matrix  mesh  respectively.  The  matrix  degrees  of  freedom  on  the  boundary 
are  eliminated  using  the  transformation. 


(72) 


< = Tqy 


(73) 


where  Tis  a transformation  matrix.  By  substituting  for  qi,"'  in  Eqn.  (72)  and  minimizing  the 
strain  energy  with  respect  to  qy  and  qi^,  we  get  the  following  relation: 
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(74) 


Figure  3 1 . Example  to  illustrate  incompatible  elements. 

(a)  problem  to  be  modeled;  (b)  matrix  mesh;  (c)  inclusion  mesh. 


Equation  (74)  can  be  solved  for  qy  and  qi^,  and  q^,'^  can  be  recovered  from  Eqn.  (73).  The 
displacement  at  any  point  in  the  unit-cell  is  obtained  by  identifying  the  inclusion  or  matrix 
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element  enclosing  the  point,  and  interpolating  the  nodal  displacements  of  that  element.  By 
repeating  this  procedure  for  several  representative  points,  the  displacement  field  within  the 
unit-cell  is  obtained. 

The  method  was  first  tested  by  predicting  the  displacement  field  for  the  unit-cell  of 
a unidirectional  composite  subject  to  uniaxial  tension.  Due  to  symmetry  of  the  unit-cell,  only 
one-quarter  of  the  unit-cell  was  modeled.  The  boundary  conditions  imposed  required  the 
edges  of  the  unit-cell  to  remain  straight  after  deformation  (Fig.  32a).  The  initial  mesh  for 
the  problem  is  shown  in  Fig.  32(b).  The  node  mismatch  at  the  inclusion-matrix  interface  can 
be  observed  in  the  figure.  Eight-node  isoparametric  elements  were  used  to  model  the  fiber 
and  matrix — 24  elements  were  used  for  the  fiber  and  98  elements  for  the  matrix.  The 
deformed  incompatible  element  configuration  is  shown  in  Fig.  32(c).  The  results  were 
compared  with  conventional  finite  element  results  by  plotting  the  displacement  fields  for 
both  cases  on  a background  25x25  mesh  (Fig.  32d).  It  was  found  that  the  incompatible 
element  displacement  field  was  identical  to  that  using  traditional  elements.  Only  one  set  of 
displacements  are  visible  in  Fig.  32(d),  because  the  displacement  fields  exactly  overlapped 
each  other. 

The  procedure  to  determine  the  background  mesh  displacement  field  is  illustrated  in 
Fig.  33.  For  a given  point  in  the  background  mesh  (Fig.  33a),  the  corresponding  incompatible 
element  containing  the  point  is  determined  (Fig.  33b).  The  displacement  for  the  mesh  point 
is  obtained  as  the  interpolation  of  nodal  displacements  of  the  found  incompatible  element 
(Fig.  33c).  A similar  procedure  is  followed  to  compute  the  background  mesh  displacements 
using  conventional  elements.  The  method  was  also  tested  by  modeling  the  transverse  shear 
of  a unidirectional  composite  (Fig.  34).  In  this  case,  one-half  of  the  unit-cell  was  modeled 
and  the  displacements  are  shown  for  a background  20x20  mesh.  In  both  examples,  it  was 
found  that  the  displacement  fields  using  incompatible  elements  and  conventional  elements 
were  identical.  A discontinuity  in  strains  in  the  vicinity  of  the  fiber-matrix  interface  was 
observed  in  both  the  examples. 
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(c)  (d) 


Figure  32.  Unidirectional  composite  subjected  to  uniaxial  tension. 

(a)  boundary  conditions;  (b)  initial  mesh  using  incompatible  elements; 
(c)  deformed  mesh;  (e)  deformed  background  mesh. 
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• background  mesh  node 
O incompatible  element  node 


Figure  33.  Calculating  displacement  field  for  the  background  mesh. 

(a)  background  mesh  showing  node  whose  displacements  are  to  be 
calculated;  (b)  incompatible  element  mesh;  (c)  matrix  element 
containing  the  background  mesh  node. 
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(c)  (d) 


Figure  34.  Unidirectional  composite  subjected  to  transverse  shear. 

(a)  boundary  conditions;  (b)  initial  mesh  using  incompatible  elements; 
(c)  deformed  mesh;  (e)  deformed  background  mesh. 
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5.1.2  Inhomogeneous  Elements  (Averaged  Gaussian  Integration) 

By  inhomogeneous  elements,  we  mean  finite  elements  which  are  comprised  of  more 
than  one  material.  Inhomogeneous  elements  in  micromechanical  analyses  were  studied  in 
detail  by  Foye  (1993).  The  advantage  of  using  inhomogeneous  elements  is  that  mesh 
generation  is  very  simple.  For  instance  a rectangular  domain  can  be  discretized  into  uniform 
rectangular  or  triangular  elements,  without  taking  into  consideration  the  constituent  material 
geometry.  The  stiffness  matrix  of  inhomogeneous  elements  represent  smeared  properties  of 
the  constituent  phases  determined  by  the  numerical  integration  scheme.  Thus  the  solution 
will  be  approximate  in  the  interface  regions.  Foye  has  developed  a modified  method  to 
evaluate  the  stiffness  matrix  of  inhomogeneous  elements  called  Replacement  Elements, 
which  predict  better  stresses  than  conventional  inhomogeneous  elements.  The  following 
details  the  procedure  to  use  inhomogeneous  elements  for  micromechanical  analysis. 

The  unit-cell  is  divided  into  uniform  rectangular  hexahedral  elements  as  shown  in 
Fig.  35.  In  general  these  elements  will  be  inhomogeneous.  The  stiffness  matrix  {K^}  for  an 
inhomogeneous  element  is  formulated  as: 

= B^CB  dV^ 

J ye 

+1  +1  +1 

B^CB  171  d^dt]d^  (75) 

-1  -1  -1 

N N N 

= X Z X 51^ 

i=\j=lk=\ 

where  is  the  domain  of  the  element,  B is  the  strain-displacement  transformation  matrix, 
C is  the  elasticity  matrix,  N is  the  number  of  Gauss  points  used  for  integration,  W is  the 
Gaussian  integration  weight  factor  and  171  is  the  determinant  of  the  Jacobian.  The  material 
property  at  each  Gauss  point  is  determined,  and  the  corresponding  elasticity  matrix 
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is  used  to  perform  the  volume  integration.  The  element  stiffness  matrix  will  thus  represent 
the  averaged  properties  of  the  constituent  materials  in  that  element. 


Figure  35.  Discretizing  a domain  using  inhomogeneous  finite  elements. 

(a)  problem  to  be  modeled;  (b)  inhomogeneous  finite  element  mesh; 
(c)  inhomogeneous  finite  element. 


85 


The  approach  was  implemented  by  writing  a code  to  identify  the  material  property 
of  an  arbitrary  Gauss  point.  The  inclusion  geometry  was  defined  by  a coarse  finite  element 
mesh.  If  the  Gauss  point  fell  within  any  of  the  inclusion  elements — the  material  property  of 
the  point  was  identified  as  that  of  inclusion.  Otherwise  the  material  property  of  the  point  was 
identified  as  matrix.  The  algorithm  to  determine  the  element  belonging  to  the  Gauss  point 
is  explained  in  the  following  section.  The  results  obtained  by  using  the  approach  to  the 
uniaxial  tension  of  a unidirectional  composite  are  shown  in  Fig.  36(a).  The  displacements 
were  in  agreement  with  those  obtained  using  homogeneous  elements  (Fig.  36b). 

The  problem  with  inhomogeneous  elements  is  that  they  cannot  represent  the  jump 
in  strains  that  can  occur  at  the  yarn-matrix  interface.  In  fact  there  are  three  strain  components 
that  can  be  discontinuous  at  the  interface,  but  the  corresponding  stresses  must  be  continuous. 
Such  a behavior  cannot  be  represented  by  inhomogeneous  elements  which  assume  a 
continuous  strain  fields  within  the  element.  This  problem  can  be  resolved  by  decomposing 
the  displacement  field  into  two  parts:  a displacement  field  qj  that  produces  a strain  field 
continuous  everywhere  in  the  unit-cell,  and  the  second  one  q2  that  has  a strain  discontinuity 
at  the  yarn-matrix  interface.  The  field  ® can  be  assumed  to  be  such  that  the  displacements 
are  identically  equal  to  zero  everywhere  in  the  matrix  and  at  the  interface,  and  exist  only  in 
the  interior  of  the  inclusions.  Thus  one  can  use  inhomogeneous  elements  for  solving  the  first 
set  of  displacements.  The  second  set  of  displacements  exists  only  in  the  inclusions,  and  they 
can  be  solved  by  discretizing  only  the  inclusion.  However  the  issue  is  determining  the 
decomposition  q=qj  +q2  ■ The  condition  for  the  decomposition  is  that  the  jump  in  interfacial 
stresses  should  be  equal  and  opposite  in  the  two  problems,  since  the  interfacial  stresses  are 
continuous  in  the  given  problem. 
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Figure  36.  Displacement  field  for  uniaxial  tension  of  a unidirectional  composite, 
(a)  inhomogeneous  elements;  (b)  homogeneous  elements. 
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Point  location  algorithm.  In  the  context  of  a textile  composite  problem,  the  material 
property  code  of  a point  is  obtained  as  follows.  The  yarn  (inclusion)  volume  is  discretized 
into  eight-node  hexahedral  isoparametric  finite  elements  as  shown  in  Fig.  37.  The  finite 
elements  are  shaped  as  triangular  prisms  (a  face  of  the  hexahedron  is  collapsed  to  an  edge). 
An  iterative  algorithm  is  used  to  determine  whether  the  point  (whose  material  property  is 
to  be  determined)  is  contained  within  the  volume  of  a given  element.  If  not,  the  procedure 
is  repeated  for  all  the  elements  in  the  yarn.  If  the  point  does  not  belong  to  any  of  the  yams, 
it  is  designated  as  a matrix  point.  If  the  point  belongs  to  a yarn  element,  the  local  yam 
direction  is  computed  as  the  direction  of  the  line  along  the  length  of  the  prismatic  element. 


Figure  37.  Finite  element  mesh  for  a yarn. 


The  algorithm  to  determine  if  a given  point  is  contained  within  the  volume  of  a finite 
element  is  explained  below.  The  procedure  is  explained  for  the  two-dimensional  case,  and 
can  be  easily  extended  to  a three-dimensional  problem. 
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(b) 
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Figure  38.  Finite  element  in: 

(a)  .xy-coordinate  system  (b)  natural  coordinate  system. 


Let  the  point,  whose  material  property  is  to  be  determined,  be  P with  coordinates  (jc,y).  It 
is  required  to  determine  whether  the  point  P is  contained  within  the  area  of  the  eight-node 
isoparametric  element  as  shown  in  Fig.  38(a).  The  element  is  mapped  into  its  natural 
coordinate  system  (Fig.  38b).  Similarly,  point  P is  mapped  to  P*in  the  element  natural 
coordinate  system.  The  element  in  the  natural  coordinate  system  is  a square  of  length 
two  units  with  edges  parallel  to  the  ^77-axes,  such  that  the  center  of  the  square  (9*  coincides 
with  the  natural  coordinate  origin.  This  means  that  we  need  to  find  whether  the  point  P*is 
contained  within  the  given  square.  This  is  possible  if  we  know  the  natural  coordinates  of 
point  P*.  However  there  does  not  exist  a transformation  from  the  jcy-coordinate  system  to 
the  natural  coordinate  system.  The  reverse  transformation  exists  as  given  below: 
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i=\ 

8 (76) 

y{^,v)  = yi 

1=1 


where  Ni  denotes  the  shape  function,  and  (xi.yi)  represents  the  jcy-coordinates  for  the  node. 

Therefore  we  use  an  iterative  procedure  to  determine  the  natural  coordinates  of  P*. 

We  substitute  ^=0  and  t/=0  in  Eqn.  (76),  and  determine  the  xy-coordinates  of  point 
O (corresponding  to  point  O*  in  the  natural  coordinate  system).  The  point  O will  not,  in 
general,  coincide  with  P.  Then  the  error,  i.e.,  the  difference  in  the  x-  and  y-coordinates 
between  points  P and  O is  calculated  as  Zlx  and  Ay  respectively.  Also  since  the  x-  and 
y-coordinates  are  functions  of  the  natural  coordinates,  Ax  and  Ay  may  be  expressed  as: 


(77) 


Equation  (77)  may  be  rewritten  in  matrix  form  as  follows: 
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where  [7(^,?;)]  is  the  Jacobian  at  (^,^).  The  Jacobian  is  computed  using  the  following 
equation: 
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Knowing  the  Jacobian,  andzl?/  can  be  evaluated  from  Eqn.  (79)  as: 
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= tel  (81) 

The  parameters  andzl?/  are  the  corrections  to  be  applied  to  the  coordinates  of  point  O* 
to  determine  the  coordinates  of  point  P*.  As  a check,  the  natural  coordinates  of  point  P*is 
then  transformed  to  its  xj-coordinates  using  Eqn.  (76),  and  compared  with  the  given 
coordinates  of  point  P.  The  procedure  is  repeated  if  the  difference  in  coordinates  is  above 
a prescribed  tolerance. 
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5.1.3  Periodic  Boundary  Conditions 


The  unit-cell  is  the  smallest  volume  element  in  the  composite  which  is  representative 
of  the  yam  architecture.  The  yam  architecture  in  the  composite  is  generally  periodic,  and  the 
unit-cell  contains  one  repeat  of  the  yarn  pattern.  Thus  the  composite  stmcture  can  be  formed 
by  assembling  the  unit-cell  in  all  three  dimensions.  When  the  composite  is  subject  to  any 
arbitrary  load,  the  microstresses  and  displacement  gradients  will  be  continuous  across  the 
faces  of  the  unit-cell.  Continuity  of  microstresses  requires  that  the  tractions  be  equal  and 
opposite  at  corresponding  points  on  opposite  faces  of  the  unit-cell.  Also  the  displacements 
between  corresponding  nodes  on  opposite  faces  of  the  unit-cell  will  differ  only  by  a constant. 


Multipoint  constraint  elements.  The  periodic  boundary  conditions  (traction  and 
displacement  boundary  conditions)  can  be  implemented  using  multipoint  constraint 
elements,  which  are  based  on  the  principle  of  Lagrange  multipliers.  The  finite  element 
formulation  with  multipoint  constraint  elements  will  be  of  the  form; 


K C 
0 


(82) 


where  {2}  is  the  matrix  of  Lagrange  multipliers  (Cook  et  al.,  1989).  However  this  method 
requires  large  storage  for  the  stiffness  matrix,  due  to  the  degrees  of  freedom  contributed  by 
the  constraint  elements  (2,).  Also  the  resulting  stiffness  matrix  will  not  be  positive  definite, 
and  the  matrix  bandwidth  will  be  very  large. 
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Transformation  equations.  The  following  section  explains  an  alternative  method  to 
enforce  the  periodic  boundary  conditions.  The  method  is  explained  for  the  two-dimensional 
case,  and  can  easily  be  extended  for  the  three-dimensional  problem.  Consider  a finite 
element  mesh  of  a unit-cell  (Fig.  39a),  where  the  periodic  boundary  conditions  are  to  be 
implemented  between  the  left  and  right  edges,  and  the  top  and  bottom  edges.  Assume  that 
each  node  has  only  one  degree  of  freedom  (d.o.f).  The  nodal  degrees  of  freedom  on  the  left 
and  bottom  edges  (except  the  prescribed  d.o.f’s)  are  designated  as  ’’master”  degrees  of 
freedom  , and  those  on  the  right  and  top  edges  as  ’’slave”  degrees  of  freedom  . The  nodal 

degrees  of  freedom  in  the  mesh  interior  are  also  classified  as  master  degrees  of  freedom. 
There  may  arise  a situation  where  the  d.o.f  for  the  corner  node  at  x/=0,  X2  = 0 is  not 
prescribed.  In  that  case,  it  is  designated  as  a master  node  with  two  slave  nodes,  whose 
coordinates  are  given  by  (x;=0,  X2  =1-2)  and  {x]=Lj,X2  = 0).  Thus  a master  d.o.f  may  have 
zero,  one  or  two  slave  d.o.f’s  depending  on  its  coordinates. 

The  slave  d.o.f’s  are  transformed  to  master  d.o.f’s  by  the  transformation: 
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qs  = Aq^  + Aqs  (83) 

where  A is  a transformation  matrix  such  that  A^J  = and  <5  denotes  the  Kronecker  delta. 

”c”  is  an  integer  function  such  that  for  the  slave  d.o.f,  c(i)  is  the  corresponding  master 
d.o.f.  Equation  (83)  can  be  rewritten  as, 

j 

— (84) 
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(A^^)j  is  the  displacement  difference  to  be  imposed  between  the  slave  d.o.f  (^^),-  and  its 
corresponding  master  d.o.f,  (qndcO)-  strain  energy  in  the  unit-cell  is  evaluated  as. 
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where  A is  the  matrix  of  Lagrange  multipliers;  Rg  and  represent  the  external  force  vectors 
for  the  slave  and  master  d.o.f ’s.  By  minimizing  the  strain  energy  with  respect  to  the  variables 
qm,  qs  and  X we  get: 
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Eliminating  qs  and  X from  Eqn.  (86),  we  arrive  at  the  expression  ; 

K*qm  = R* 


(86) 


(87) 


where 


K*  = [Kmm  + A^KsA  + A'^Ksm  + Krr,A] 
R = Rm  + A^Rs  — [Kms  + A^Kss]Aqs 


(88) 
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Equation  (87)  can  be  solved  for  q,n  i recovered  from  Eqn.  (83).  Thus  the  number  of 
unknowns  in  the  linear  equation  solver  are  reduced,  and  also  the  resulting  stiffness  matrix 
(.S’*)  is  positive-definite. 


5.2  Mesh  Generation 

Currently  approximate  methods  such  as  averaged  Gaussian  integration.  Selective 
Averaging  Method  and  Isostrain  Method  are  popular  because  of  their  ease  of  computation 
and  requirement  of  minimum  computer  storage  and  time.  However  the  advancements  in 
computer  hardware  and  computational  technology  will  make  it  possible  to  use  a large 
number  of  degrees  of  freedom  for  micromechanical  analysis.  Then  there  will  be  a need  for 
mesh  generation  techniques  for  creating  homogeneous  elements  in  a unit-cell.  In  the 
following  we  describe  two  methods  for  discretizing  the  unit-cell  with  homogeneous 
elements. 

5.2.1  Node  Migration  Method 

The  unit-cell  is  meshed  with  tetrahedral  solid  elements  in  an  arbitrary  fashion.  This 
will  be  called  the  primitive  mesh.  The  elements  are  identified  as  homogeneous  or 
inhomogeneous.  Then  the  nodes  of  the  inhomogeneous  elements  are  moved  to  the  interface 
by  using  an  heuristic  algorithm.  In  each  inhomogeneous  element,  the  node  closest  to  the 
interface  is  allowed  to  migrate  to  the  interface  first.  After  each  cycle,  some  elements  will 
transform  into  homogeneous  elements.  Then  the  process  is  repeated  until  all  the  elements 
become  homogeneous.  The  mesh  thus  generated  is  called  the  intermediate  mesh. 

The  intermediate  mesh  will  have  some  elements  distorted  due  to  node  migration. 
This  distortion  can  be  removed  by  subjecting  the  mesh  to  an  annealing  process,  by  which 
the  distortions  concentrated  near  the  interfaces  are  distributed  among  other  elements  also. 
An  example  of  this  method  is  depicted  in  Fig.  40.  The  purpose  here  is  to  discretize  a square 
unit-cell  containing  two  circular  inclusions.  The  figure  shows  the  initial  mesh  containing 


94 


uniform  square  elements,  the  intermediate  mesh  of  homogeneous  elements  and  the  annealed 
mesh. 
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5.2.2  Modified  Node  Migration  Method 

In  this  method  an  initial  mesh  is  created  as  in  the  previous  method.  The  elements  are 
divided  into  three  groups:  homogeneous  matrix  elements,  homogeneous  yam  elements  and 
inhomogeneous  elements  which  contain  both  matrix  and  yarn.  Then  the  inhomogeneous 
elements  are  redesignated  as  matrix  elements,  thus  leaving  only  two  kinds  of  elements.  We 
will  call  this  as  the  intermediate  mesh.  Thus  the  finite  element  representation  of  the  yarn  will 
be  smaller  than  the  actual  yarn,  i.e.,  the  yarn  mesh  will  be  fully  contained  within  the  actual 
yam.  Then  the  nodes  on  the  yarn-matrix  interface  in  the  intermediate  mesh  are  allowed  to 
migrate  to  the  nearest  point  on  the  actual  yarn-matrix  interface.  As  before  a 
three-dimensional  finite  element  program  with  fictitious  material  properties  is  used  in  this 
step  to  obtain  an  annealed  mesh.  An  example  of  the  final  mesh  for  a square  unit-cell  with 
two  inclusions  in  the  shape  of  a quarter-circle  is  shown  in  Fig.  41. 


CHAPTER  6 

CONCLUSIONS  AND  SUGGESTIONS  FOR  FUTURE  WORK 


The  previous  chapters  presented  some  of  the  aspects  involved  in  the  modeling  of 
textile  structural  composites.  Micromechanical  analyses  were  developed  to  predict  the 
three-dimensional  elastic  constants  and  CTE’s  for  textile  composite  materials.  The  effect  of 
stress  gradients  in  thin  textile  composites  was  highlighted,  and  consequently,  an  independent 
micromechanical  analysis  was  developed  for  such  composites.  The  issues  involved  in 
modeling  a thin  textile  composite  were  demonstrated  by  first  modeling  the  composite  as  a 
beam  and  then  as  a plate.  Two  codes  called  ^TE%-\Q  and  fiTEx-20  were  developed  in 
FORTRAN  77,  which  implemented  the  finite  element  procedure  presented  in  Chapter  2,  and 
the  SAM  procedure  presented  in  Chapter  3 respectively.  The  predicted  macroscale 
thermo-elastic  constants  compared  very  well  with  existing  models  for  textile  composite 
materials. 

The  finite  element  analysis  and  SAM  procedures  assume  that  the  unit-cell  of  the 
composite  is  shaped  as  a rectangular  hexahedron.  Both  analyses  can  be  easily  modified  for 
other  unit-cell  geometries.  For  example,  in  the  finite  element  procedure,  the  relations  for  the 
periodic  boundary  conditions  will  change  based  on  the  unit-cell  geometry.  However  the 
assumption,  that  the  unit-cell  is  rectangular,  limited  the  textile  composite  examples  for 
which  the  above  codes  can  be  implemented.  For  instance,  the  unit-cell  for  a braided 
composite  is  hexahedral  with  the  included  angle  between  the  edges  equal  to  the  braid  angle. 
The  braid  angle,  in  general,  will  not  be  equal  to  a right  angle.  Thus  the  effective  composite 
properties  could  be  computed  only  for  woven  architectures,  where  the  unit-cell  is  a 
rectangular  hexahedron. 
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A considerable  amount  of  time  was  invested  in  developing  mesh  generation 
algorithms  and  alternative  finite  element  procedures,  in  the  context  of  modeling  textile 
composites.  Some  of  these  ideas  are  presented  in  Chapter  5.  The  mesh  generation  algorithms 
such  as  the  node  migration  method  were  very  effective  for  two-dimensional  examples. 
However  for  the  three-dimensional  case,  some  of  the  elements  near  the  interface  collapsed 
to  create  degenerate  and  skewed  elements.  Further,  extremely  small  tetrahedral  elements 
were  required  to  capture  the  interstitial  matrix  geometry.  Computer  storage  and  CPU  time 
limitations  made  the  use  of  such  a finite  element  mesh  impractical.  Initially  multipoint 
constraint  elements  were  used  to  impose  periodic  boundary  conditions  in  ^TEx-10.  To 
reduce  computer  storage,  the  transformation  method  (for  periodic  boundary  conditions)  and 
the  skyline  solver  (Bathe,  1982)  were  incorporated  in  the  finite  element  code.  Further 
reduction  in  computer  storage  can  be  realized  by  using  sparse  matrix  solvers  instead  of  the 
skyline  solver. 

The  averaged  Gaussian  integration  technique  was  used  to  compute  the 
inhomogeneous  element  stiffness  matrix  infiTEx-lO.  The  inhomogeneous  element  averages 
the  properties  of  the  constituent  materials  in  the  element,  and  assumes  that  the  strain 
distribution  is  continuous  within  the  element.  Consequently,  the  inhomogeneous  element 
cannot  capture  the  discontinuity  in  strains  across  the  yarn-matrix  interface.  Thus  the 
corresponding  microstresses  computed  in  the  vicinity  of  the  interface  will  be  inaccurate.  The 
accuracy  of  the  microstresses  may  be  improved  by  refining  the  (inhomogeneous  element) 
mesh,  but  that  would  greatly  increase  the  degrees  of  freedom  in  the  numerical  model.  This 
emphasizes  the  need  for  effective  finite  element  preprocessing  codes  for  meshing  the 
unit-cell  of  the  textile  composite  with  homogenous  elements.  It  will  also  be  useful  to 
interface  the  finite  element  results  with  suitable  graphics  software.  This  will  enable  the  user 
to  get  a feel  for  the  problem  by  directly  visualizing  the  unit-cell  deformed  configurations, 
and  the  microstress/microstrain  distributions  within  the  unit-cell. 
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Chapter  4 presented  failure  envelope  results  for  a plain- weave  textile  beam.  Implicit 
to  the  procedure  for  determining  the  failure  envelope,  was  the  assumption  that  microscopic 
failure  in  the  unit-cell  translated  to  macroscopic  failure  of  the  composite.  The  procedure  is 
simplistic  and  conservative,  but  yet  of  use  to  a structural  designer.  More  work  remains  to  be 
done  in  obtaining  results  for  failure  envelopes  using  the  strength  models  for  textile 
composite  continuum  and  for  a composite  plate,  as  explained  in  Chapter  4.  An  issue  which 
this  study,  due  to  time  limitations,  does  not  address  is  parametric  analyses  to  study  the  effect 
of  changing  constituent  material  properties  and  fiber  volume  fraction  on  composite 
properties.  Both  the  finite  element  analysis  and  SAM  can  be  easily  extended  for  other 
macroscopic  composite  properties  such  as  thermal  conductivities,  electro-magnetic 
properties  and  so  on. 
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